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phase  of  the  frequency  response  of  mixed  phase  signals  and  systems. 

A  number  of  techniques  are  applied  to  estimation  of  the  phase  fre¬ 
quency  response  of  the  speech  production  tract  from  quasi -periodic 
speech  segments.  Methods  of  phase  estimation  are  categorized  as 
indirect  or  direct.  A  subject  of  the  indirect  procedures  yield  a 
closed  form  solution  for  retrieving  the  phase  from  the  magnitude 
of  a  mixed  phase  frequency  response  and  a  priori  knowledge  about 
the  corresponding  signal.  Linear  iterative  algorithms  are  also 
developed  for  retrieving  the  phase  from  the  magnitude,  and,  simi¬ 
larly,  the  magnitude  from  the  phase,  with  a  causality  or  finite 
duration  constraint  imposed  on  the  desired  signal.  The  iterative 
algorithm  for  magnitude  retrieval  provides  an  alternative  to  the 
Hilbert  transform  for  obtaining  the  magnitude  from  the  phase  of  a 
minimum  phase  signal,  but  without  the  need  of  an  unwrapped  phase. 

In  addition,  it  serves  as  the  major  component  within  a  new  phase 
unwrapping  algorithm  which  does  not  require  modulo  2^  considerations 
An  alternate  indirect  strategy  changes  a  phase  estimation  problem 
to  one  of  magnitude  estimation  by  modifying  a  quasi-periodic  wave¬ 
form  so  that  the  desired  impulse  response  takes  on  a  minimum  phase 
characteri sti c Di rect  approaches  rely  on  harmonic  samples  of  a 
frequency  response,  or  of  the  principal  value  of  its  phase.  Specif¬ 
ically,  time-domain  windowing  and  frequency-domain  interpolation 
are  applied  to  a  quasi-periodic  waveform  in  estimating  the  unwrappe  I 
phase  at  harmonics.  Mixed  phase  estimates  from  these  direct  ap¬ 
proaches  are  incorporated  within  homomorphic  and  spectral  envelope 
speech  analysis-synthesis  systems  with  a  second-order  improvement 
in  quality  over  their  minimum  phase  counterparts. 
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ABSTRACT 

This  dissertation  addresses  the  problem  of  estimating  the  phase  of 
the  frequency  response  of  mixed  phase  signals  and  systems.  A  number  of 
techniques  are  applied  to  estimation  of  the  phase  of  the  frequency  response 
of  the  speech  production  tract  from  quasi -periodic  speech  segments.  Methods 
of  phase  estimation  are  categorized  as  indirect  or  direct.  A  subset  of  the 
Indirect  procedures  yield  a  closed  form  solution  for  retrieving  the  phase 
from  the  magnitude  of  a  mixed  phase  frequency  response  and 'a  priori  Ncnow- 
ledge  about  the  corresponding  signal.  Linear  Iterative  algorithms  are 
also  developed  for  retrieving  the  phase  from  the  magnitude,  and,  similarly, 
the  magnitude  from  the  phase,  with  a  causality  or  finite  duration  con¬ 
straint  Imposed  on  the  desired  signal.  The  Iterative  algorithm  for  magni¬ 
tude  retrieval  provides  an  alternative  to  the  Hilbert  transform  for  obtain¬ 
ing  the  magnitude  from  the  phase  of  a  minimum  phase  signal,  but  without  the 
need  of  an  unwrapped  phase.  In  addition.  It  serves  as  the  major  component 
within  a  new  phase  unwrapping  algorithm  which  does  not  require  modulo  Zv 
considerations.  An  alternate  Indirect  strategy  changes  a  phase  estimation 
problem  to  one  of  magnitude  estimation  by  modifying  a  quasi -periodic  wave¬ 
form  so  that  the  desired  Impulse  response  takes  on  a  minimum  phase  charac¬ 
teristic.  Direct  approaches  rely  on  harmonic  samples  of  a  frequency  re¬ 
sponse,  or  of  the  principal  value  of  Its  phase.  Specifically,  time-domain 
windowing  and  frequency-domain  Interpolation  are  applied  to  a  quasi -periodic 
waveform  In  estimating  the  unwrapped  phase  at  harmonics.  Mixed  phase  esti¬ 
mates  from  these  direct  approaches  are  incorporated  within  homomorphic  and 
spectral  envelope  speech  analysis-synthesis  systems  with  a  second-order 
Improvement  In  quality  over  their  minimum  phase  counterparts. 
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CHAPTER  1 
INTRODUCTION 

In  many  physical  situations,  we  encounter  the  problem  of  recovering  a 
signal  which  Is  not  directly  accessible  to  measurement,  but  for  which  par¬ 
tial  knowledge  of  Its  Fourier  transform  can  be  determined.  This  partial 
knowledge  may,  under  certain  conditions,  be  sufficient  to  recover  the 
entire  Fourier  transform.  For  example,  when  a  signal  Is  known  a  priori 
to  be  minimum  or  maximum  phase,  the  Fourier  transform  can  be  recovered 
within  a  sign  or  scale  factor  from  only  Its  magnitude  or  phase[22]. 

Within  this  dissertation,  we  shall.  In  particular,  consider  techniques 
which  use  various  kinds  of  partial  Information  In  estimating  the  phase  of 
the  Fourier  transform  of  mixed  phase  discrete-time  signals  and  systems. 

A  number  of  these  techniques  are  applied  specifically  to  the  analysis  of 
voiced  speech.  That  Is,  we  consider  estimation  of  the  phase  of  the  fre¬ 
quency  response  of  a  signal  which  consists  of  the  convolution  of  the  com¬ 
bined  vocal  and  nasal  tract  impulse  response,  and  the  glottal  wavelet[5]. 

A  simplified  model  of  voiced  speech  sounds  such  as  vowels,  nasalized  vowels, 
and  nasal  consonants  is  given  by  the  convolution  of  this  response  with  a 
train  of  equally  spaced  Impulses  representing  the  periodicity  of  the  vocal 
cord  excitation  funct1on[5].  This  type  of  representation,  l.e.  a  "quasl- 
perlodlc"  waveform, Is  also  characteristic  of  other  waveforms  found  for 
example  In  biology,  music,  and  many  other  acoustic  disturbances. 

We  shall  apply  a  subset  of  our  procedures  for  phase  estimation  to  two 
speech  analysis -synthes is  systems:  (1)  the  homomorphic  system  proposed  by 
0ppenhe1m[23],  and  (11)  the  spectral  envelope  system  proposed  by  Paul[26]. 
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The  original  schemes  Introduced  either  a  zero  or  minimal  phase  Impulse 
response  derived  from  a  magnitude  estimate  only. 

Although  many  of  the  phase  estimation  procedures  we  consider  are 
developed  In  the  context  of  speech  analysis,  our  techniques  have  more 
general  Implications  and  potential  use  In  other  areas  where  phase  estima¬ 
tion  Is  required.  For  example,  a  number  of  results  may  be  extended  to 
two-dimensional  signals  and  systems,  and  thus  are  possibly  applicable  to 
such  areas  as  Image  or  seismic  signal  processing. 

1.1  Scope  of  Thesis 

Most  deconvolution  techniques  for  recovering  an  Impulse  response  from 
a  quasi -periodic  waveform  are  directed  primarily  to  estimating  the  magni¬ 
tude  of  the  desired  frequency  response.  Linear  pred1ct1on[2,15]  and  homo¬ 
morphic  f11ter1ng[23,24],  for  example,  have  been  applied  quite  successfully 
to  magnitude  estimation. 

The  phase  estimation  procedures  we  shall  consider  are  roughly  classi¬ 
fied  as  either  "Indirect"  or  "direct".  Our  class  of  indirect  schemes 
capitalizes  on  either  knowing  the  desired  magnitude  or  being  capable  of 
deriving  an  accurate  estimate  of  it  by  conventional  deconvolution  tech¬ 
niques.  Clearly,  when  the  desired  signal  is  minimum  or  maximum  phase,  the 
phase  can  be  obtained  by  applying  a  Hilbert  transform  to  the  logarithm  of 
the  given  magnitude.  More  generally,  we  shall  use  various  kinds  of  a 
priori  Information  about  the  desired  signal  along  with  the  magnitude  to 
unambiguously  determine  a  mixed  phase  function.  In  particular,  imposing 
causality  or  a  finite  length  constraint  on  the  signal  and  specifying  a  few 
samples  of  the  phase,  or  the  first  few  points  of  a  discrete-time  sequence, 
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In  some  cases.  Is  sufficient  to  uniquely  characterize  the  entire  phase 
function.  For  example,  as  we  shall  see,  for  a  certain  class  of  causal 
sequences  only  the  Initial  value  of  the  sequence  Is  necessary. 

An  alternate  Indirect  means  of  recovering  phase  from  magnitude  con¬ 
verts  the  phase  estimation  problem  to  a  magnitude  estimation  problem  by 
modifying  the  speech  waveform  so  that  the  desired  Impulse  response  takes 
on  a  minimum  phase  characteristic.  Specifically,  a  class  of  Invertible 
transformations  Is  derived  which  are  suitable  to  changing  the  general  prob¬ 
lem  of  deconvolution.  Involving  both  magnitude  and  phase  estimation  of  a 
mixed  phase  sequence  from  a  quasl-perlodlc  waveform,  to  a  deconvolution 
problem  where  only  a  magnitude  estimate  Is  required.  This  procedure  Is 
applied  to  homomorphic  deconvolut1on[24].  The  sensitivities  Inherent  In 
this  deconvolution  scheme,  due  to  the  requirement  of  an  unwrapped  phase 
[33],  are  therefore  avoided.  In  addition.  It  appears  that  this  indirect 
approach  to  phase  estimation  by  homomorphic  deconvolution  is  less  sensitive 
to  noise  disturbances  than  a  direct  approach  which  requires  an  unwrapped 
phase.  Since  most  noise  reduction  systems  estimate  only  the  magnitude  of 
a  frequency  response,  this  technique  Is  also  potentially  applicable  to 
signal  enhancement [14],  where  both  magnitude  and  phase  estimates  are  ob¬ 
tained  through  a  magnitude  estimate  only. 

A  linear  Iterative  algorithm  Is  also  developed  for  retrieving  the 
phase  from  the  magnitude  and  a  priori  Information  about  a  desired  signal. 
The  algorithm  obtains  the  phase  by  Iteratively  Imposing  the  known  magnitude 
function  In  the  frequency  domain,  and  a  priori  Information  about  the  signal 
in  the  time  domain.  The  algorithm  therefore  falls  within  a  class  which 
encompasses,  for  example,  the  recently  proposed  iterative  techniques  by 


n 

Gerchberg[6]  and  Papoul1s[2S]  for  extrapolation  of  a  finite  portion  of  a 
bandllmlted  signal. 

An  analogous  Iterative  procedure  Is  likewise  presented  for  retrieving 
the  magnitude  of  the  frequency  response  from  the  phase.  When  the  sequence 
is  minimum  phase,  applying  the  Hilbert  transform  to  the  unwrapped  phase  re¬ 
covers  the  logarithm  of  the  magnitude  (within  a  scale  factor)[22].  Our 
Iterative  algorithm,  however,  when  Imposing  the  a  priori  knowledge  of  caus¬ 
ality,  recovers  the  magnitude  from  only  the  principal  value  of  the  phase  of 
a  minimum  phase  sequence.  The  procedure  of  phase  unwrapping  Is  therefore 
avoided.  The  algorithm  under  this  particular  constraint  also  serves  as  the 
major  component  within  a  new  phase  unwrapping  algorithm  which  does  not  re¬ 
quire  the  conventional  modulo  2it  considerations [32],  As  we  shall  see,  an 
unwrapped  phase  Is  needed  not  only  In  phase  estimation  by  homomorphic  decon¬ 
volution,  but  also  In  other  techniques  within  the  thesis. 

In  parallel  with  this  thesis  development,  Hayes  et  al[10]  have  demon¬ 
strated  that  when  Imposing  a  finite  length  constraint  on  a  sequence,  the 
magnitude  of  Its  frequency  response  In  general  Is  uniquely  specified  by 
the  phase  (within  a  scale  factor).  As  a  result,  we  show  that  our  Iterative 
algorithm  can,  as  well,  recover  the  magnitude  from  the  phase  of  a  finite 
length  mixed  phase  sequence. 

The  direct  approaches  we  shall  consider  do  not  require  an  estimate 
of  a  magnitude  function,  but  rather  require  partial  knowledge  of  Its  phase 
function  which  Is  derived  from  the  phase  of  the  speech  waveform.  Homomor¬ 
phic  deconvolution  Is  one  such  approach. 

An  alternate  direct  strategy  addresses  the  problem  of  phase  estima¬ 
tion  from  partial  knowledge  consisting  of  harmonic  samples  of  the  desired 
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frequency  response  or  of  Its  principal  phase  value.  Our  estimation  algor¬ 
ithms  attempt  to  preserve  the  unwrapped  phase  at  harmonics  —  a  procedure 
referred  to  as  "phase  tracking".  Phase  tracking  ensures  that  the  un¬ 
wrapped  phase  estimate  at  harmonics  equals  the  unwrapped  phase  of  the  de¬ 
sired  frequency  response.  Such  a  property  Is  desirable  since  under  a  band- 
limited  constraint  on  the  unwrapped  phase,  the  entire  unwrapped  phase  may 
be  approximately  recovered  from  only  the  given  harmonic  samples.  One 
method  of  phase  tracking  Invokes  an  Interpolation  procedure  on  harmonic 
samples,  or  samples  of  the  principal  phase  value.  In  particular,  we  derive 
conditions  on  the  desired  unwrapped  phase  under  which  a  simple  linear  Inter¬ 
polation  scheme  across  two  consecutive  samples  preserves  the  unwrapped 
phase  at  harmonics. 

This  study  leads  to  a  heuristic  understanding  of  the  Interaction 
between  windowing  a  speech  waveform  modeled  as  exactly  periodic  over  a 
short  time,  and  the  nature  of  the  unwrapped  phase  of  the  windowed  waveform. 
In  particular,  constraints  on  the  duration  and  position  of  a  specific  class 
of  windows  are  derived  for  guaranteeing  phase  tracking  by  the  windowing 
procedure  itself. 

Phase  estimates  derived  from  the  techniques  of  phase  tracking  are  in¬ 
corporated  within  our  two  speech  analysis-synthesis  schemes.  Constraints 
on  time-domain  windowing  play  a  major  role  In  governing  the  accuracy  of 
the  phase  estimate  within  these  systems.  In  the  homomorphic  scheme  tailor¬ 
ing  the  duration  and  position  of  the  window  Is  used  to  Improve  the  phase 
estimate  by  homomorphic  filtering. 

Linear  Interpolation  In  the  frequency  domain  can  likewise  be  viewed 
In  the  time  doamln  as  multiplication  by  a  window.  The  position  of  this 
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window  Is  also  Important  in  obtaining  an  accurate  unwrapped  phase  estimate 
at  harmonics.  This  consideration  will  be  made  In  designing  a  high  quality 
spectral  envelope  speech  analysis-synthesis  system  with  a  phase  estimate 
derived  from  our  linear  Interpolation  procedure. 

With  an  alternative  Interpretation  of  windowing  both  analysis-synthesis 
schemes  are  shown  to  be  Identity  systems  with  respect  to  reconstruction  of 
a  periodic  waveform. 

Informal  listening  tests  Indicate  a  small,  but  perceptible  Improvement 
In  "quality"  when  In  these  systems  a  mixed  phase  reconstruction  replaces 
Its  minimum  phase  counterpart. 

1.2  Phase  In  Speech 

The  phase  of  the  frequency  response  of  the  speech  production  tract  has 
generally  been  considered  less  Important  than  the  magnitude  function  In 
generating  high  quality  synthetic  speech  within  speech  analysis-synthesis 
systems.  Experiments,  however,  have  been  reported  demonstrating  that  the 
envelope  of  a  periodic  waveform  can  be  an  Important  factor  In  audible  per- 
cept1on[17].  In  particular,  it  Influences  sensations  of  roughness  or 
smoothness.  The  change  from  roughness  to  smoothness  may  be  accomplished 
by  changing  the  phase  or  magnitude  of  a  particular  component  or  set  of  com¬ 
ponents  of  a  periodic  waveform,  and  the  degree  of  roughness  Is  related  to 
the  relative  length  and  depth  of  the  recurrent  depressions  In  the  envelope. 
This  "peaklness"  factor  Is  thus  a  determinant  In  roughness  of  a  periodic 
waveform. 

The  ultimate  goal  In  any  speech  analysis-synthesis  system  Is  to  ex¬ 


tract  from  the  speech  waveform  In  analysis  perceptually  Important  informa¬ 
tion  which  Is  used  in  synthesis  to  reconstruct  the  original  waveform. 
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Short-tin*  s pooch  segments  or*  loosoly  categorized  os  olthor  volcod  or 
unvoiced.  Tho  Input-output  system  model  for  volcod  segments  wos  glvon  In 
tho  Introduction.  Unvolcod  s pooch  Is  llkowlso  modeled  ovor  o  short-time 
os  tho  convolution  of  the  combined  vocol  ond  nosol  troct  Impulse  response 
with  o  whit*  nolso  oxcltotlon.  The  onolysls  stop*  usuolly  extracts  tho 
magnitude  of  tho  Fourier  transform  of  tho  doslrod  Impulse  response,  o 
volcod/ unvolcod  decision,  and  a  pitch  measurement  for  voiced  signals. 

Cither  a  zero  phase  or  a  minimum  phase  function  Is  1ntroduced[2,23]. 

Both  the  zero  and  the  minimum  phase  Impulse  response  estimates  are 
characterized  by  peaklness.  The  minimum  phase  estimate,  for  example, 
yields  a  maximum  energy  concentration  at  the  signal's  origin  due  to  Its 
minimum  delay  property[22].  Although  there  does  not  exist  a  demonstration 
of  a  correlation  between  peaklness  and  quality  degradation  In  speech 
analysis-synthesis  systems,  the  highly  quasi -periodic  characteristic  of 
speech  and  perceptual  tests  on  the  envelope  of  a  periodic  waveform  suggests 
the  possibility  of  "quality"  Improvement  by  reduction  of  peaklness  with  an 
accurate  phase  estimate  of  the  desired  response.  Two  results  developed  In 
parallel  with  this  thesis  support  this  conjecture.  First,  Atal  and  David 
[1]  have  shown  that  the  quality  of  Linear  Prediction  Coding  (LPC)  Is  Im¬ 
proved  by  Introducing  an  approximation  to  the  LPC  error  residual.  For 
voiced  speech  the  prediction  residual  Is  quasl-perlodlc  with  the  same  per¬ 
iod  as  the  speech  waveform.  A  pitch  period  long  segment  of  the  prediction 
residual  can  be  expressed  In  terms  of  a  Fourier  series  expansion  as  a  sum 
of  contributions  of  the  fundamental  and  the  Individual  harmonics.  The 
contribution  of  a  particular  harmonic  (e.g.  the  kth  harmonic)  Is  given  by 
a  cosine  function  of  a  particular  frequency,  amplitude,  and  phase: 
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A^  co*(iukn  +  ek).  The  Fourier  series  representation  allows  the  variation 
of  the  amplitude  and  phase  of  different  harmonics  In  any  desired  fashion. 
Listening  tests  Indicate  a  second  order,  but  perceptible  change  In  speech 
quality  for  any  magnitude  condition  when  the  phase  changes  from  zero  to 
the  phase  derived  at  each  harmonic.  Furthermore,  there  Is  only  a  slight 
difference  between  the  original  phase  and  a  fixed  frequency-dependent  phase. 
This  fixed  phase  was  generated  by  computing  the  medium  group  delay  at  each 
harmonic  over  all  pitch  periods  Included  In  a  sentence-like  structure. 

Go1d[8]  In  preliminary  listening  tests  has  also  found  quality  Improve¬ 
ment  by  introducing  a  phase  function  within  the  channel  vocoder.  A  fixed 
frequency-dependent  mixed  phase  was  Introduced  along  with  the  phase  function 
of  a  minimum  phase  vocal  tract  Impulse  response  estimate  derived  from  LPC. 

1.3  Outline  of  Chapters 

The  thesis  begins  In  Chapter  2  with  developing  a  framework  for  the 
phase  estimation  problem  In  the  context  of  speech  analysis  and  discrete¬ 
time  signals  and  systems.  A  review  Is  given  of  a  specific  set  of  results 
to  be  used  throughout  the  following  chapters. 

In  Chapter  3  we  consider  the  Indirect  approach  of  estimating  the  phase 
from  the  magnitude  of  a  Fourier  transform.  Constraints  and  solutions  are 
presented  for  unambiguous  retrieval  of  the  phase  from  the  magnitude  func¬ 
tion.  These  constraints  are  then  applied  within  a  linear  Iterative  algor¬ 
ithm  for  phase  retrieval. 

In  the  second  portion  of  this  chapter,  we  Investigate  methods  of  con¬ 
verting  the  phase  estimation  problem  to  a  magnitude  estimation  problem. 

A  specific  class  of  Invertible  transformations  Is  derived  for  changing  an 
arbitrary  mixed  phase  sequence  to  a  minimum  phase  sequence,  and  which  are 
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suitable  to  preserving  the  convolutional  characteristic  of  a  quasi -periodic 
waveform. 

In  Chapter  4  we  first  review  homomorphic  deconvolution  for  directly 
estimating  a  phase  function  from  a  quasl-perlodlc  waveform.  We  next  de¬ 
scribe  the  heuristics  of  the  sensitivity  of  this  approach  due  to  the  need 
of  an  unwrapped  phase  function.  Homomorphic  deconvolution  without  an  un¬ 
wrapped  phase  Is  then  used  as  an  Illustration  of  the  Indirect  approach  by 
transformation  of  a  quasl-perlodlc  waveform  given  In  Chapter  3.  Finally 
In  this  chapter,  we  present  a  new  phase  unwrapping  algorithm  which  does  not 
require  modulo  2ir  considerations.  Toward  this  objective,  we  first  discuss 
a  class  of  constraints  which  guarantees  unambiguous  magnitude  retrieval 
(within  a  scale  factor)  from  the  phase  of  a  Fourier  transform.  A  linear 
Iterative  algorithm,  the  dual  to  that  In  Chapter  3,  Is  then  developed  to 
retrieve  the  magnitude  from  the  phase  under  such  constraints.  We  Illustrate 
the  use  of  this  algorithm  by  recovering  the  magnitude  of  a  minimum  phase 
sequence  (which  yields  the  same  result  as  would  a  Hilbert  transform),  and 
by  recovering  also  the  magnitude  of  a  finite  length  mixed  phase  sequence. 
Finally,  this  Iterative  technique  Is  used  as  the  major  component  In  our  new 
phase  unwrapping  algorithm. 

In  Chapter  5  we  take  a  direct  approach  to  phase  estimation  within  a 
different  context  than  found  In  Chapters  3  and  4.  Specifically,  we  assume 
only  harmonic  samples  of  the  desired  frequency  response  are  available. 
Techniques  of  phase  tracking  are  developed,  which  Involve  both  linear  Inter¬ 
polation  In  the  frequency  domain,  and  windowing  In  the  time  domain. 

In  Chapter  6  techniques  of  phase  tracking  derived  In  Chapter  5  are 
applied  In  designing  the  two  previously  described  high  quality  speech 


analysis-synthesis  systems  with  mixed  phase.  Finally  In  this  chapter 
results  of  Informal  listening  tests  are  presented,  where  mixed  and  minimum 
phase  reconstructions  are  compared. 

Lastly,  in  Chapter  7  a  summary  of  the  main  results  of  the  thesis  Is 
given.  Me  also  suggest  a  direction  and  some  potential  areas  of  future 
research. 
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CHAPTER  2 

A  FRAMEWORK  FOR  PHASE  ESTIMATION  IN  SPEECH  ANALYSIS 

In  this  chapter  a  framework  Is  developed  for  the  phase  estimation 
problem  In  speech  analysis  In  the  context  of  discrete-time  signals  and 
systems.  We  first  present  the  general  problem  of  estimating  the  magnitude 
and  phase  of  the  Fourier  transform  of  a  sequence  h(n)  which  Is  convolved 
with  a  train  of  equally  spaced  pulses  p(n):  x(n)  »  h(n)  *  p(n).  This  for¬ 
mulation  Is  applicable  to  voiced  speech  segments  modeled  over  a  short  dur¬ 
ation  by  a  "quasi -periodic"  waveform  which  Is  produced  by  exciting  the 
vocal  and  nasal  tract  with  pulses  of  air,  l.e.,  the  glottal  wavelet,  caused 
by  vibration  of  the  vocal  cords[5].  We  shall  assume  that  the  desired  se¬ 
quence  h(n)  consists  of  the  convolution  of  the  Impulse  response  of  the  vocal 
and  nasal  tract  and  glottal  wavelet. 

The  phase  estimation  component  of  this  deconvolution  problem  Is  often 
tied  to  magnitude  estimation.  In  fact,  for  a  certain  class  of  sequences 
the  phase  function  can  be  derived  from  the  magnitude  which  Is  often  more 
easily  and  directly  measurable  than  phase.  We  are  therefore  led  to  a  review 
of  minimum  and  maximum  phase  sequences  which  fall  within  this  class,  and 
for  which  the  Hilbert  transform  relations  exist  between  the  continuous 
phase  and  logmagnltude  of  the  Fourier  transform.  These  relations  enable 
the  phase  to  be  reconstructed  from  the  magnitude  and,  likewise,  the  magni¬ 
tude  from  the  phase.  In  general,  without  additional  constraints  neither 
the  magnitude  nor  the  phase  Is  sufficient  to  completely  characterize  a 
mixed  phase  sequence. 
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Since  the  phase  of  a  complex  number  Is  a  multivalued  function.  It  Is 
generally  not  unique.  With  appropriate  constraints,  however,  a  unique  func¬ 
tion  can  be  derived  and  Is  termed  the  unwrapped  phase.  Besides  being  a  re¬ 
quirement  In  the  Hilbert  transform  relatlonsfor  minimum  and  maximum  phase 
sequences,  the  unwrapped  phase  Is  also  an  Important  component  In  the  theory 
of  homomorphic  deconvolution,  and  In  other  direct  phase  estimation  procedures. 

Our  primary  objective  In  this  chapter  Is  to  collect  a  specific  set  of 
results  which  are  useful  In  the  chapters  to  follow.  We  will  avoid  a  detailed 
discussion  of  subtle  mathematical  Issues,  and  present  only  those  results 
necessary  In  formulating  certain  techniques  of  phase  estimation. 

2.1  Deconvolution  of  a  Quasi -Periodic  Waveform 

Consider  a  sequence  x(n)  which  consists  of  the  convolution  of  a  sequence 
h(n)  and  p(n)  a  train  of  equally  spaced  pulses: 

x(n)  *  h(n)  *  p(n)  (2.1) 

where  p(n)  Is  given  by 

p(n)  *  l  a(k)  5(n-kP)  (2.2) 

k 

and  where  P  Is  the  pulse  spacing  (l.e. ,  the  "pitch  period"  In  the  context 
of  voiced  speech),  <5(n)  Is  the  unit-sample  sequence[22],  and  k  may  range 
over  finite  or  Infinite  extent. 

From  physical  considerations  h(n)  Is  assumed  to  be  causal,  l.e.,  h(n)»0 

for  n<0,  and  stable,  l.e.,  £  |h(n)|  <  •  .  The  z-transform  of  h(n)  Is  given 

n 

H(z)  ■  l  ti(n)z'"  (2.3) 

n*0 

and  has  a  region  of  convergence  which  encompasses  the  unit  circle  and  the 
entire  z-plane  outside  the  unit  circle.  Including  z  ■  •  .  H(z)  Is  referred 
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to  as  the  system  function,  and  when  evaluated  on  the  unit  circle  yields  the 
frequency  response  of  the  system:  H(z)(  2«eXp[ju]*H(u)*  In  general,  H M  Is 
complex  and  can  be  expressed  In  terms  of  Its  real  and  Imaginary  parts  as 

H(w)  -  Hr(«)  +  j  Ht (ui)  (2.4) 

where  r  and  1  denote  real  and  Imagin'  respectively.  In  terms  of  Its  mag¬ 
nitude  and  phase  (i.e. ,  polar  form),  (2.4)  is  expressed  as 

H(a>)  -  |H(u) |  exp[jeh(a>)]  (2.5) 

where  eh(o»)  *  arg[H(u)]. 

With  X(z),  H(z),  and  P(z)  denoting  the  z- transforms  of  x(n),  h(n),  and 
p(n),  respectively.  It  follows  from  (2.1)  that  X(z)  can  be  written  as 

X(z)  *  H(z)  P(z)  (2.6) 

X(z)  represents  the  z-transform  of  a  stable  sequence  when  p(n)  is  absolutely 

sumnable,  I.e.,  £  |p(n)|  <  •  ,  and  in  particular  when  p(n)  Is  bounded  and 
n 

of  finite  extent  so  that  the  range  of  k  Is  finite. 

The  Fourier  transform  of  x(n)  is  found  by  evaluating  X(z)  on  the  unit 
circle  and  in  polar  form  is  expressed  by 

X(uj)  *  X^l^expi;^]  a  |X(«)|  exp[j6x(o.)] 

*  |H(«)  |  |P(«)  |  exp[jeh(u»)]  exp[j0p(tu)]  (2.7) 

where,  exU)  *  arg[X(u»)]  and  ep(u>)  *  arg[0p(w)]. 

The  goal  of  deconvolution  of  the  quasi -periodic  sequence  of  (2.1)  Is 
to  extract  the  magnitude  and  phase  of  K(u)  and  P(u>).  Estimation  of  |H(u)| 
has  been  extensively  Investigated  and  Is  generally  an  easier  problem  than 


21 
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estimation  of  eh(w).  Extraction  of  P(<d)  Is  often  not  required  In  speech 
analysis  but  rather  only  the  period  P  is  sought  —  a  problem  denoted  as 
"pitch  detection".  Within  this  dissertation,  we  shall  restrict  ourselves  to 
recovering  the  phase  of  H(w). 

2.2  Minimum.  Maximum,  and  Mixed  Phase  Sequences 

In  this  section,  we  formulate  the  definitions  and  properties  of  minimum, 

maximum,  and  mixed  phase  sequences.  The  unwrapped  phase  function  Is  also 

defined  and  its  characteristics  are  discussed. 

2.2.1  Definitions 

The  definitions  of  minimum  and  maximum  phase  sequences  are  completely 
analogous  to  each  other  and  are  given  below. 

Minimum  Phase  Sequences:  A  complex  function  H(z)  of  a  complex  variable 
z  is  minimum  phase  if  it  is  analytic  and  its  reciprocal  H-1(z)  is  also 

analytic  for  |z|  ^  1  in  the  z-plane.  A  minimum  phase  sequence  is  then  de¬ 

fined  as  a  sequence  whose  z- transform  is  minimum  phase.  It  follows  from 
this  definition  that  a  necessary,  but  not  sufficient,  condition  for  h(n)  to 
be  minimum  phase  is  that  it  be  causal,  stable,  and  nonzero  at  n  =  0.  That 
is, 

h(n)  3  0  n  <  0  (2.8.a) 

h(n)  f  0  n  »  0  (2.8.b) 

l  |h(n) |  <  »  (2.8.c) 

n 

Maximum  Phase  Sequences:  A  complex  function  H(z)  of  a  complex  variable 
z  is  maximum  phase  if  it  is  analytic  and  its  reciprocal  H'^z)  is  also  analy¬ 
tic  for  | z [  <.1  In  the  z-plane.  A  maximum  phase  sequence  is  defined  as  a 
sequence  whose  z-transform  is  maximum  phase.  Such  sequences  are  anti-causal. 
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stable,  and  nonzero  at  n  *  0.  That  is. 


h(n)  *  0 

n  >  0 

(2.9. a) 

h(n)  *  0 

n  *  0 

(2.9. b) 

I  |h(n) |  <  » 

(2.9. c) 

n 


Conditions  2. 9. a,  b,  and  c,  however,  are  necessary,  but  not  sufficient,  for 
a  maximum  phase  characteristic. 

Mixed  Phase  Sequences:  A  complex  function  H(z)  of  a  complex  variable 
z  which  is  neither  minimum  nor  maximum  phase  is  termed  mixed  phase.  A  mixed 
phase  sequence  has  a  z-transform  which  is  mixed  phase.  Is  stable,  and  may 
or  may  not  be  causal  or  anti -causal. 

2.2.2  Hilbert  Transform  Relations  for  Minimum  Phase  Sequences 

A  consequence  of  the  properties  of  a  minimum  phase  sequence  is  a  useful 
relationship  between  the  logarithm  of  |H(w)|  and  eh (cu) .  Specifically,  the 
Hilbert  transform  relations [22]  provide  a  means  of  retrieving  the  phase  from 
the  magnitude,  and  the  magnitude  (within  a  scale  factor)  from  the  phase  for 
such  sequences. 

To  derive  these  relations  consider  the  complex  logarithm  of  H(z)  given 
by 

H(z)  *  log[H(z)] 

3  log|H(z) |  +  j  arg[H(z)]  (2.10) 

If  H(z)  is  viewed  as  the  z-transform  of  a  real  sequence  h(n),  then  when  h(n) 
Is  causal,  h(n)  can  be  completely  recovered  from  its  even  component 
he(n)  *  (h(n)  +  h(-n))/2,  or  its  odd  component  hQ(n)  *  (h(n)  -  h(-n))/2  and 
h(0)[22].  That  is,  h(n)  is  given  by 
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or 


(V“> 

n  a  0 

/ . 

2He(n) 

n  >  0 

to 

n  <  0 

fi(0> 

n  «  0 

n  >  0 

to 

n  <  0 

(2.11a) 


(2.11b) 


The  Fourier  transform  of  hg(n)  is  log  |H(<u)|  ,  and  the  Fourier  transform  of 

A 

hQ(n)  is  arg[H(w)].  Therefore,  a  consequence  of  (2.11a)  and  (2.11b)  is  that 
the  imaginary  component  of  H(w),  arg[H(<a)]  can  be  recovered  from  its  real 
component,  log  |H(a>)|.  Likewise,  the  real  component,  log | H(o>)|  can  be  re¬ 
covered  from  its  imaginary  component,  arg[H(a>)]  within  the  additive  con¬ 
stant  h(o).  In  fact,  it  is  possible  to  obtain  direct  relations  between 
log|H(u>)|  and  arg[H(u>)],  i.e.,  the  Hilbert  transform  relations.  Neverthe¬ 
less,  magnitude  and  phase  reconstruction  is  generally  easier  to  perform 
indirectly  through  (2.11a)  and  (2.11b). 

The  requirement  that  log|H(w)|  and  arg[H(u>)]  be  a  Hilbert  transform 
pair  Is  often  referred  to  as  the  minimum  phase  condition.  It  corresponds 
to  the  requirement  that  the  sequence  h(n)  be  causal  and  stable  and  can  be 

shown  to  be  equivalent  to  the  definition  of  a  minimum  phase  sequence  given 

* 

in  the  preceding  section.  The  analogous  case  with  h(n)  anti -causal  and 
stable  corresponds  to  the  maximum  phase  condition. 


i 
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2.2.3  The  Unwrapped  Phase  Function 

A  A 

For  h(n)  of  the  preceding  section  to  be  causal  and  stable,  H(z)  must  be 
analytic  In  the  region  |z|  >.1.  In  considering  this  analytlcltyy  we  must 
appropriately  define  arg[H(z)]  since  any  multiple  of  2ir  can  be  added  to  the 
phase  without  affecting  the  value  of  H(z),  and  thus  arg[H(z)]  in  general 
will  be  discontinuous. 

The  phase  is  therefore  ambiguous  to  within  a  2ir  multiple,  and  can  be 
expressed  by 

arg[H(z)]  *  ARG[H(z)]  +  2irq  (2.12a) 

where  q  ■  0,1 ,2,. . .  and 

-ir  <  ARG[H(z)]  <  ir  (2.12b) 

ARG[H(z)]  Is  termed  the  principal  value  of  arg[H(z)].  This  ambiguity  Is 

A 

resolved  by  the  fact  that  analyticity  of  H(z)  implies  that  Its  real  and 
imaginary  parts  must  be  continuous  functions  of  z,  and  consequently  if 

A 

H(z)  Is  to  be  analytic,  we  must  define  arg[H(z)]  to  be  a  continuous  func¬ 
tion.  Furthermore,  since  h(n)  is  assumed  real,  arg[H(z)]  will  be  defined 
so  that  for  z  a  exp[ju]  it  is  odd,  periodic  in  u>  with  period  2ir,  and  a  con¬ 
tinuous  function  of  w[32].  A  phase  function  which  satisfies  these  proper¬ 
ties  is  termed  the  unwrapped  phase  function. 

One  approach  to  computing  the  unwrapped  phase  Is  to  assume  that  a 
continuous  phase  function  Is  obtained  by  Integration  of  the  phase  deriva¬ 
tive.  Assuming  a  differentiable  complex  logarithm  in  (2.10),  and  evaluating 
the  logarithmic  derivative  on  the  unit  circle,  we  obtain 
H'(u)  -  H'  (ui)/H(u) 

-  h;u)  +  j  Hju) 


(2.13) 
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where  the  prime  denotes  differentiation  with  respect  to  «. 
obtain 

Hj(w)  ■  arg'[H(w)] 

Hr(«)H;.(o>)  - 
Hf2(«)  + 


From  (2.13)  we 


(2.14) 


Integrating  (2.14)  ensures  that  arg[H(u>)]  is  a  continuous  function  of 
u.  In  addition,  to  ensure  that  arg[H(u)]  be  odd  and  periodic  In  u  with 
period  2ir,  the  following  two  conditions  must  be  met: 


Since, 


and 


arg[H(ui)] 


O) 


0,1T 


H(oj) 


Hi  *  0 


l  h(n) 
n 


arg[H(u)] 


u 


IT 


ir 

arg'[H(cu)]  dw 

. 

0 


(2.15) 


(2.16) 


(2.17) 


It  follows  that  only  sequences  with  a  positive  mean  and  a  zero  mean  phase 
derivative  are  compatible  with  the  above  requirements. 

n 

We  define  the  linear  phase  component  of  H(z)  as  z  ,  where 
nQ  *  arg[H(ir)]/ir  and  where  arg[H(ir)]  Is  given  In  (2.17).  Thus,  it  is  per¬ 
haps  necessary  that  h(n)  be  Inverted  and  shifted  (removing  the  linear  phase 
component  corresponds  to  a  shift  [22])  In  order  that  (2.15)  be  satisfied. 

Nevertheless,  it  will  be  useful  to  modify  our  definition  of  unwrapped 

n 

phase  to  encompass  a  linear  phase  component.  With  the  presence  of  z  ,  the 
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unwrapped  phase  Is  defined  over  the  half-open  Interval  [o.ir)  by  Integration 
of  the  phase  derivative,  and  Is  continuous  In  this  Interval.  However,  to 
be  odd.  It  must  be  discontinuous  at  ir[32]. 

Since  many  of  the  techniques  developed  In  this  thesis  will  be  Implement¬ 
ed  on  a  digital  computer,  discrete  Fourier  transform  (DFT)  Implementations 
are  required.  In  particular,  we  need  to  determine  samples  of  the  continuous 
phase  function.  One  technique  of  obtaining  samples  of  the  unwrapped  phase 
first  computes  the  principal  value  of  the  unwrapped  phase  using  Inverse 
tangent  routines  and  then  "unwraps"  by  simply  adding  appropriate  multiples 
of  2ir  to  the  principal  value  until  the  discontinuities  Induced  by  the 
modulo  2ir  operation  are  removed[32].  A  more  recent  technique  has  been  pro¬ 
posed  that  combines  the  Information  In  both  the  phase  derivative  and  the 
principal  value  of  the  phase  Into  an  adaptive  numerical  Integration  unwrap¬ 
ping  scheme[34].  The  method  adds  appropriate  multiples  of  2ir  to  the  princi¬ 
pal  value  of  the  phase  until  the  phase  is  "consistent"  with  the  numerically 
integrated  phase  derivative. 

2.3  Sequences  with  Rational  z-Transforms 

In  this  section  we  restrict  ourselves  to  sequences  whose  z- transform  is 

given  by  a  rational  function  of  the  form 

m1  "o 

n  (l-akz-1)  n  (1-b.  z) 

H(z)  -  Az  0  {£* - ^ -  (2.18) 

n  (1-c. z_1 )  n  (l-dkz) 

k-1  K  k»l  K 

n 

where  |ak|,  |bk|,  |ck|,  and  | dk|  are  less  than  or  equal  to  unity  and  z  Is 
the  linear  phase  component.  Factors  of  the  form  (1-a^z’1)  and  (1-c^z-1) 
correspond  to  zeros  and  poles  on  or  Inside  the  unit  circle,  and  the  factors 


A 
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(1-b^z)  and  (1-d^z)  correspond  to  zeros  and  poles  on  or  outside  the  unit 
circle. 

2.3.1  Minimum  and  Maximum  Phase  Rational  z-Transforms 

We  shall  Interpret  the  results  of  the  previous  sections  In  the  context 

A 

of  this  more  restrictive  class  of  sequences.  Since  H(z)  of  (2.10)  (or 

equivalently  H(z)  and  H"1 (z) )  must  be  analytic  In  the  region  |z|  >,  1,  for 

H(z)  to  be  a  minimum  phase  function,  there  can  be  no  poles  or  zeros  of  H(z) 

on  or  outside  the  unit  circle.  This  includes  poles  and  zeros  at  Infinity; 

l.e. ,  If  H(z)  Is  minimum  phase  11m  H(z)  must  be  a  nonzero  finite  constant. 

z-*» 

The  Implication  of  this  restriction  Is  that  factors  of  the  form  (l-bkz)  and 
(l-dkz)  do  not  exist,  and  a  linear  phase  component  Is  not  present. 

A  In  (2.8)  may  be  positive  or  negative  when  H(z)  Is  minimum  phase. 
However,  for  H(z)  to  be  compatible  with  the  requirement  that  there  exists  an 
odd,  continuous  phase  function,  A  must  be  positive. 

The  z-transform  of  a  normalized  (i.e.,  A*l)  minimum  phase  sequence  can 


therefore  be  expressed  as 


Wz> 


n  (1-a.z  ') 
k*l  * 

Pi  -1 

n  (l-ckz 

k=1  * 


(2.19) 


A  completely  analogous  formulation  can  be  made  for  maximum  phase 
sequences  for  which  poles  and  zeros  lie  outside  the  unit  circle.  The  cor¬ 
responding  normalized  z-transform  is  given  by 


Hmax<z> 


n  (1-b.z) 
k»1  K 

P0 

n  (1-dkz) 
k*l  * 


(2.20) 
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Therefore,  a  mixed  phase  sequence  with  no  zeros  on  the  unit  circle 
Is  expressed  by 

h(n)  »  Ahmin(n)  *  hmx(n)  *  5  (n-nQ)  (2.21) 

where  hm^n(n)  and  hTOX(n)  correspond  to  (2.19)  and  (2.20),  and  «s(n-nQ) 
accounts  for  a  linear  phase  component.  A  mixed  phase  sequence  In  general, 
however,  may  have  zeros  on  the  unit  circle. 

2.3.2  Magnitude- Phase  Relations  for  Mixed  Phase  Rational  z-Transforms 

In  this  section  we  consider  the  problem  of  reconstructing  the  magni¬ 
tude  from  the  phase,  and  the  phase  from  the  magnitude  of  a  Fourier  trans¬ 
form,  for  a  mixed  phase  rational  z-transform. 

Suppose  we  are  given  the  magnitude  of  the  Fourier  transform  of  a 
causal  sequence.  A  phase  function  can  be  derived  Indirectly  through  (2.11a) 
to  obtain  a  minimum  phase  z-transform  denoted  by  Hmp(z).  The  poles 
(assumed  inside  the  unit  circle)  of  Hrap(z)  are  the  poles  of  the  original 
z-transform  H(z),  and  zeros  inside  the  unit  circle  are  also  left  Intact. 
However,  zeros  of  H(z)  outside  the  unit  circle  are  mapped  to  their  con¬ 
jugate  reciprocal  locat1ons[22].  The  original  z-transform  H(z)  can 
therefore  be  represented  by  the  cascade  of  a  minimum  phase  system  and 
an  all-pass  system,  A(z):  H(z)  »  Hmp(z)A(z),  where  an  all-pass  Is  de¬ 
fined  as  a  system  for  which  the  magnitude  of  the  frequency  response  Is 
unity  for  all  u>.  In  particular,  an  arbitrary  all-pass  A(z)  can  be  shown 
to  consist  of  a  cascade  of  factors  of  the  form: 

1  -  a*z 
1  -  az*1 

where  j  a |  <  1  . 


(2.22) 
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Consequently,  such  systems  have  the  property  that  their  poles  and  zeros 
occur  at  conjugate  reciprocal  locations. 

Since  the  original  sequence  h{n)  Is  causal,  knowledge  of  the  magni¬ 
tude  of  H(u)  uniquely  specifies  the  poles  of  H(.z)  which  from  above  are 
equal  to  the  poles  of  H  (z).  Therefore,  with  knowledge  of  the  existence 

Hip 

M 

of  M  zeros  of  H(z)  there  exists  a  maximum  of  2  different  phase  functions 
(excluding  linear  phase)  for  a  given  magnitude  function.  These  phase  func¬ 
tions  can  be  generated  by  reflecting  zeros  about  the  unit  circle  through 
(2.22).  Without  additional  a  priori  knowledge  the  exact  location  of  the 
zeros  and  thus  the  original  phase  cannot  be  determined. 

Consider  now  the  dual  problem  of  recovering  the  magnitude  from  the 
phase  of  the  z- transform  on  the  unit  circle.  We  assume  the  linear  phase 
component  is  removed  to  satisfy  the  phase  continuity  condition  for  the 
Hilbert  transform  relations.  We  shall  see  that  the  ambiguity  in  obtaining 
the  magnitude  function  Is  of  a  different  nature  from  the  problem  of  phase 
retrieval  from  magnitude. 

As  before,  the  minimum  phase  counterpart  to  the  original  z-transform 
is  found  by  applying  (2.11b).  To  see  the  resulting  pole-zero  pattern 
consider  the  all-pass  function 

Hap(z)  -  H(z)/H(z'1)  (2.23) 

Evaluating  H._(z)  on  the  unit  circle,  we  obtain 
ap 

I,  •  ,xp[ju]  ■  «PWV“>]  <2-24> 


where  eh(u)  is  the  known  phase  function.  (2.24)  can  be  written  as 


■  Wz>  • 

(2.25) 

where  z-exp[j<*»] 

*Vin^  *  Hmin^z^Hmax^z  ^ 

(2.26a) 

and. 

**max^  *  Hmax^z^Hm1n^z  ^ 

(2.26b) 

Since  Hm1n(z)  and  Hmax(z)  both  have  Identical  phase  9h(o>),  on  the 
unit  circle,  Hml.n(z)  corresponds  to  the  minimum  phase  counterpart  of  H(z) 
with  phase  0^(01).  From  (2.26a),  the  poles  and  zeros  Inside  the  unit  cir¬ 
cle  of  H(z)  remain  Intact,  while  the  zeros  outside  the  unit  circle  have 
been  reflected  as  poles  to  their  conjugate  reciprocal  locations. 

The  original  z-transform  H(z)  can  therefore  be  represented  by  the 
cascade  of  a  minimum  phase  system  and  a  zero-phase  system  B(z)  for  which 

the  phase  of  the  frequency  response  Is  zero  for  all  w:  H(z)  *  H  (z)B(z). 

min 

In  particular,  an  arbitrary  zero-phase  system  B(z)  can  be  shown  to  consist 
of  a  cascade  of  factors  of  the  form 

[(1  -  az”1 ) (1  -  a*z)]  i1  .  (2.27) 

Given  that  there  exist  no  poles  outside  the  unit  circle  and  that  the  ori¬ 
ginal  z-transform  has  a  total  of  L  poles  Inside  the  unit  circle  or  zeros 
outside  the  unit  circle,  there  exists  a  maximum  of  2**  different  magnitude 
functions  for  a  given  phase  function.  These  magnitude  functions  can  be 


generated  by  reflecting  poles  to  zeros  and  zeros  to  poles  through  the 
zero-phase  function  of  (2.27).  For  example,  a  maximum  phase  zero,  repre< 
sented  by  the  factor  (1-bz),  where  |b|  <  1  ,  Is  reflected  to  a  minimum 


phase  pole  by  the  operation: 

,  -1  ,  -1 
(1  -  b*z"‘)  •  (1  -  bz ) C ( 1  -  b*z  )(1  -  bz)] 

(2.28) 

Therefore,  without  additional  a  priori  knowledge,  the  magnitude  function 
cannot  be  uniquely  specified  from  the  phase. 


i 

i 


t 

I 
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CHAPTER  3 

PHASE  RETRIEVAL  FROM  MAGNITUDE 

We  have  seen  In  the  preceding  chapter  that  for  minimum  or  maximum 
phase  sequences  the  phase  of  a  Fourier  transform  Is  uniquely  recoverable 
from  Its  magnitude.  When  a  sequence  does  not  fall  within  this  class.  It 
Is  reasonable  to  seek  alternative  a  priori  knowledge  about  the  sequence 
which  Is  sufficient  to  unambiguously  retrieve  the  phase  from  the  magnitude 
function. 

For  example,  we  may  know  a  priori  that  the  zeros  of  a  rational  z- 
transform  H(z)  lie  outside  the  unit  circle  and  that  the  corresponding 
sequence  Is  causal  and  nonzero  at  the  origin.  The  zeros  of  the  minimum 
phase  z-transform,  with  the  given  magnitude  on  the  unit  circle,  from  (2.22), 
fall  at  conjugate  reciprocal  locations  to  those  of  the  original  z-transform. 
H(z)  and  consequently  the  phase  of  H(u>)  can  therefore  be  retrfeved[10]. 

In  this  chapter  we  shall  consider  two  categories  of  constraints: 

(I)  Constraints  on  values  of  points  of  the  sequence 

(II)  Constraints  on  the  values  of  samples  of  the  phase  function 

In  a  number  of  cases  these  constraints  lead  to  a  set  of  linear  equations 
for  obtaining  parameters  of  a  rational  z-transform,  and  thus  an  Indirect 
means  of  retrieving  the  phase  function. 

An  alternative  method  of  phase  retrieval  Invokes  a  linear  Iterative 
procedure  which  Incorporates  known  values  of  the  sequence  and  the  known  mag¬ 
nitude  function.  This  Iterative  algorithm  Is  particularly' useful  when  a 
linear  solution  does  not  exist  for  the  given  constraints. 

Recently  there  has  been  a  great  deal  of  Interest  In  Iterative  algor- 
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Ithms  for  signal  reconstruction.  Gerchberg[6]  and  Papoul1s[25],  for 
example,  have  developed  an  iteration  for  the  determination  of  a  bandlimited 
continuous-time  signal  x(t)  in  terms  of  a  finite  segment  of  x(t).  Gerchberg 
and  Saxton[7]  utilize  the  magnitude  of  a  complex  signal  and  the  magnitude 
of  its  Fourier  transform  in  an  iterative  fashion  to  obtain  the  phase  func¬ 
tions  In  both  the  time  and  frequency  domains.  Mersereau  and  Schafer[20] 
also  have  recently  performed  a  comparative  study  of  various  Iterative  decon¬ 
volution  algorithms.  These  techniques  Impose  a  time-limited,  band-limited, 
and  positivity  constraint  on  the  output  of  the  deconvolution. 

Such  algorithms  and,  likewise,  the  algorithm  proposed  within  this 
chapter  all  fall  within  a  class  where  information  about  the  sequence  and 
Its  Fourier  transform  Is  Iteratively  Imposed.  In  the  next  chapter,  we 
describe  an  Iterative  algorithm,  also  of  this  class,  which  recovers  the  mag¬ 
nitude  of  a  Fourier  transform  from  its  phase. 

When  constraints  of  the  above  kind  are  not  specified  or  accuracy  of 
the  magnitude  Is  uncertain,  we  consider  the  alternative  procedure  of  con¬ 
verting  a  phase  estimation  problem  to  a  magnitude  estimation  problem  by 
modification  of  the  speech  waveform.  The  second  main  area  of  this  chapter. 
In  particular.  Investigates  methods  to  transform  a  mixed  phase  sequence  to 
a  minimum  phase  sequence.  With  an  appropriate  Invertible  transformation  a 
magnitude  function  is  then  sufficient  to  completely  characterize  the  modi¬ 
fied  sequence.  Such  transformations  are  useful  In  changing  the  general 
problem  of  deconvolving  a  mixed  phase  sequence  from  a  quasi -periodic  wave¬ 
form  to  a  deconvolution  problem  where  only  an  estimate  of  the  magnitude  of 
the  Fourier  transform  of  the  modified  sequence  Is  required.  The  Invertlbll- 
Ity  of  the  transformation  allows  approximate  recovery  of  the  desired  phase 
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function. 

A  second  area  of  potential  application  of  this  technique  Is  In  signal 
enhancement[14].  We  can  modify  a  sequence  to  take  on  a  minimum  phase  char¬ 
acteristic  before  degradation  by  additive  noise.  Most  noise  reduction 
systems  such  as  spectral  subtraction  estimate  only  the  magnitude  of  the 
Fourier  transform[13].  Transformation  of  a  sequence  allows  both  magnitude 
and  phase  estimation  through  only  a  magnitude  estimate. 


3.1  Phase  Retrieval  From  Magnitude  with  Constraints 

When  a  sequence  h(n)  is  causal,  stable  and  has  a  rational  z-transform. 


H(z)  Is  given  by 


H(z) 


n  mi  mo 

A z°  n  (1-a.z*1)  n  (1-b.z) 
k»l  K  k»l  * 


r1  i 

“  (l-ckz  ') 
k*l 


(3.1) 


where  |a^|,  |b|J,  |CjJ,  and  |djJ  are  less  than  unity.  We  argued  in  section 

no 

2.3.2  that  (neglecting  the  linear  phase  component  z  )  there  exists  a  maxi- 
mum  of  2  possible  phase  functions  when  |H(w)|  is  given,  where  M»m^+mQ  is 
the  number  of  zeros  In  (3.1).  We  shall  now  consider  constraints  on  h(n)  and 
e^U)  to  resolve  this  phase  ambiguity.  These  constraints  are  not  exhaust¬ 
ive,  but  Indicative  of  the  nature  of  requirements  for  phase  recovery. 

3.1.1  Constraints  on  the  Sequence 

One  method  of  guaranteeing  a  unique  phase  function  is  to  constrain 
values  of  points  of  the  sequence.  In  section  3. 1.1.1,  we  specify  the  first 
M+l  points  of  h(n)  to  remove  phase  ambiguity.  Alternatively,  In  section 

3. 1.1.2  we  demonstrate  that  for  a  restricted  class  of  sequences,  the  Ini- 
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tial  value  of  h(n),  h(0)  Is  sufficient  for  a  unique  phase  determination. 

3. 1.1.1  Infinite  Length  Sequences  and  the  Pade  Approximation 

The  first  method  we  consider  borrows  the  philosophy  of  the  Pade 
approximation  of  a  rational  z-tranform[18].  In  this  technique  the  para¬ 
meters  of  the  rational  model  are  chosen  to  exactly  match  the  first 
p^  +  m^  +  mQ  +  1  points  of  the  sequence,  where  p^ ,  m^,  and  mQ  are  given 
in  (3.1).  In  a  similar  manner,  we  shall  show  that  when  the  first  m^+mQ+l 
points  of  the  sequence  are  matched  and  the  magnitude  is  given,  the  para¬ 
meters  of  the  rational  model  can  be  uniquely  determined,  and  thus  a  phase 
function  unambiguously  specified. 

Assuming  that  h(0)^0,  and  expanding  the  numerator  and  denominator 
functions  of  (3.1)  as  polynomials  in  z-1,  we  obtain 

H(z)  -  N(z)/D(z) 

M  N 

=  l  a(k)z'k/  l  b(k)z’k  (3.2) 

k=0  k=0 

where  M*m^+mo  and  N=p...  When  the  magnitude  of  H(u>)  is  known,  from  section 
2.3.2,  the  denominator  polynomial  D(z)  of  H(z)  can  be  determined [18].  Cross 
multiplying  by  D(z),  we  can  express  (3.2)  in  the  time  domain  as  a  convolu¬ 
tion: 

a(n)  =  b(n)  *  h(n) 

N 

=  l  b(k)h(n-k)  (3.3) 

k=0 

A 

Given  the  coefficients  b(n)  and  h(n)  for  n*0,l ,2,. . .M,  the  numerator  coef¬ 
ficients  a(n)  are  easily  computed  from  (3.3). 
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3.1. 1.2  Knowledge  of  h(0) 

When  h(n)  is  causal  and  h(0)^0,  H(z)  In  (3.1)  is  expressed  by 


"M  -1  "'o  .1 

A  n1  O-a.z  ')  n  (z  -bJ 


Pi  1 
n  (l-ckz  J 


(3.4) 


The  zeros  of  (3.4)  may  be  flipped  inside  and  outside  the  unit  circle  by  an 
all -pass  function  A(z)  consisting  of  factors  represented  by 


[(z"1  -  a*)/(l  -  a  z "Mr 


(3.5) 


This  particular  form  maintains  |H(ai)|  ,  while  preserving  the  causality  of 
the  original  sequence  h(n).  (3.5)  also  implies  a  nonzero  initial  value. 

Therefore,  from  the  Initial  Value  Theorem[22],  for  each  zero  geometry, 
with  | H (oj) |  held  constant,  the  initial  value  h(OV  is  given  by 

h (0 )  3  lim  H(z)A(z) 


A  n  (-aj  n  (-b.) 
kel  K  keO  * 


(3.6) 


where  A(z)  is  an  all-pass  function,  and  where  I  represents  the  set  of  in¬ 
tegers  which  correspond  to  maximum  phase  factors  of  the  form  (z'^-a^), 
and  0  represents  the  set  of  Integers  which  correspond  to  maximum  phase 
factors  of  the  form  (z‘^-bk).  Thus  given  arbitrary  sets  I-j  and  0^,  and 
I2  and  0^,  each  possible  zero  geometry  corresponds  to  a  distinct  h(0) 
under  the  constraint^!  ]  n  (_a  )  n  (-b  ) 


i  n  (-aj  n  (-b.)  (3.7) 

kel2  K  ke02  K 
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The  constraints  represented  by  (3.7)  and  knowledge  of~hT0)~h«vfr  neither  a 
simple  meaningful  interpretation,  nor  a  linear  solution,  as  does  matching 
the  first  M+l  points  of  h(n).  Nevertheless,  they  provide  a  flavor  for  the 
requirements  necessary  for  phase  retrieval. 

In  summary,  when  (3.7)  is  satisfied,  |H(<u)|  is  given,  and  H(z)  is  of 
the  form  in  (3.4),  a  unique  phase  function  is  guaranteed  for  each  distinct 
h(0). 

3.1.2  Constraints  on  The  Phase  Function 

Consider  the  situation  where  the  magnitude  is  known  and  samples  of  the 
phase  function  are  given.  We  shall  show  that  when  the  number  of  samples  of 
the  phase  in  the  frequency  interval  [0,ir]  is  greater  than  the  order  of  the 
numerator  M  in  (3.1),  the  phase  is  uniquely  specified. 

Given  samples  of  arg(H(z))  on  the  unit  circle  at  ,.  ..u^, 

samples  of  H(z)  can  also  be  determined  on  the  unit  circle  since  (H(uj)|  is 
known.  Cross-multiplying  by  D(z)  in  (3.2)  with  z  *  exp[juim],  we  obtain 
M 

l  a(k)exp[-jku)m]  =  H(u)m)D(u>m)  (3.8) 

k=0 

With  M+l  samples  of  the  phase  function  (3.8)  represents  M+l  equations  in 
M+l  unknowns.  The  real  component  of  this  set  of  complex  linear  equations 
is  given  in  matrix  form  by 
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or, 

Ca  -  Nr  (3.9b) 

where  Nr(m)  «  Re(H(oim)D(a>m)). 

McClellan  and  Parks[19]  have  demonstrated  that  the  set  of  functions  (l,cos(w), 
cos(2u)...cos(Mw)}  is  a  Chebyshev  set  for  u>e[0,Ti].  An  implication  of  this 
result  Is  that  when  0  <  u>_  <  it  and  the  samples  u>  are  distinct,  the  matrix 
C  is  invertible.  There  exists,  therefore,  a  unique  set  of  numerator  coef¬ 
ficients  a(k)  and  consequently  an  unambiguous  phase  function. 

In  summary,  when  |H(ui)|  is  known,  h(n)  Is  causal,  and  M*1  distinct 
phase  samples  are  given  In  the  interval  [0,ir],  the  phase  function  is  unique. 
Clearly,  a  consequence  of  this  result  is  that  any  finite  segment  of  the 
phase  function  is  also  sufficient  to  remove  ambiguity  since  M+l  phase 
samples  can  always  be  obtained  from  such  a  segment. 

3.2  An  Iterative  Procedure  to  Retrieve  Phase  From  Magnitude 

In  this  section  a  linear  Iterative  procedure  is  developed  for  phase 
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retrieval.  The  method  Invokes  values  of  the  sequence  and  the  magnitude 
of  the  Fourier  transform  which  guarantee  an  unambiguous  phase  function 
when  the  z-transform  Is  assumed  rational. 

3.2.1  Theory 

The  Iterative  algorithm  will  now  be  described  and  a  proof  Is  given 
showing  that  a  defined  error  must  decrease  monotonlcally  as  the  algorithm 
Iterates. 

The  algorithm  Is  Illustrated  in  Figure  3.1.  We  begin  with  an  Initial 
guess  9q(o>)  of  the  desired  phase,  and  Inverse  transform  the  function 
M(a>)exp[j9g(u))]t  where  M(u))*|H(aj) |  Is  the  given  magnitude.  This  step  yields 
hg(n),  the  initial  estimate  of  h(n).  Next,  the  known  values  of  h(n)  for 
nel  (where  I  is  the  set  of  Integers  for  which  h(n)  Is  given)  are  Incorpor¬ 
ated  in  the  Initial  estlmte  hQ(n)  to  obtain  hg(n).  The  magnitude  of  the 
Fourier  transform  of  hg(n)  Is  then  replaced  by  the  given  magnitude  and  the 
procedure  is  repeated.  The  steps  involved  in  one  iteration  are  sunmarlzed 
below.  hk(n),  9k(ui),  and  Mk(u)  are  the  signal,  phase,  and  magnitude  esti¬ 
mates,  respectively  on  the  kth  iteration  and  hk+^(n)  is  defined  by 


hk<"> 

nil 

h(n) 

nel 

(3.10) 


I 


Fig.  3.1  Iterative  algorithm  to  recover  phase 
from  magnitude 


Iteration  to  Recover  Phase  From  Magnitude 

(i)  Inverse  transform  M(w)exp[jek(w)]:hk(n) 

(ii)  Replace  hk(n)  with  h(n)  for  nel:fik+^(n) 

(iii)  Forward  transform  hk+-j(n):Mk+^(&»)exp[j0k+1(a))] 

(iv)  Replace  Mk+1(w)  by  M(w):M(u)exp[j9k+1  (o>)] 

(v)  Repeat 

The  above  steps  complete  one  iteration  and  the  algorithm  Is  continued 
for  as  many  desired  number  of  Iterations.  In  order  to  demonstrate  that  the 
algorithm  results  in  a  "reasonable"  phase  estimate,  we  wish  to  choose  an 
error  function  which  Is  monotone  decreasing  through  each  pass  of  the  algor¬ 
ithm.  One  error  function  we  have  considered  Is  the  sum  of  the  squared  dif¬ 
ferences  between  h(n)  and  the  estimate  hk(n)  on  each  iteration: 

Ek  *  I  |h(n)  -  hk(n)|2  (3.11) 

Clearly,  from  (3.10)  the  error  decreases  in  the  time  domain  whenever 
hk(n)^h(n)  for  nel,  and  stays  the  same  whenever  hk(n)=h(n)  for  nel. 

That  is, 

Ek  >  I  |h(n)  -  hk+1(n)|2  .  (3.12) 

n 

However,  In  the  frequency  domain  it  is  not  possible  to  show  through  a 
vector  argument  that  |H(w)  -  Hk+1(u>)|  monotonically  decreases  at  each  fre¬ 
quency  cd  when  the  known  magnitude  is  Incorporated  In  Hk+i(u>),  the  Fourier 
transform  of  hk+-j  (n ) . 

For  this  reason  we  choose  an  alternative  error  function  which  can  be 
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easily  shown  to  monotonlcally  decrease  upon  each  pass  of  the  algorithm.  The 
error  function  Is  defined  as  the  mean  squared  difference  between  the  known 
magnitude  and  the  estimate  Mk(u)  on  each  Iteration: 


Ek'?7 


|M(u)  -  M^(u) |  du 


-ir 


We  shall  show  In  two  steps  that  Ek  Is  monotone  decreasing: 


(3.13) 

Ek  -  ^+1  • 


Error  Reduction  In  the  Time  Domain 

With  the  Identity  |exp[jek(u)]|2  *  1,  the  expression  for  Ek  In  (3.13) 
can  be  written  as 


E 


k 


1 

27 


■IT 

|M(oj)  -  Mk(u)|2  | exp[j9k (a»)] |2  du 

j 

-ir 


ir 

|M(«)  exp[jek(u)]  -  Mk(u)  exp[jek(u)]|2  du 

-IT 


|Hk(«)  -  Hk(u)|2  du  (3.14) 

where  Hk(u)  and  Hk(u)  are  the  Fourier  transforms  of  hk(n)  and  hk (n ) ,  res¬ 
pectively.  From  Parseval's  Theorem  (3,14)  Is  given  In  the  time  domain  by 


Ek  ”  I  lhk(n>  “  \(">T 

n 


From  (3.10)  It  follows  that 


hk(n)  -  hk+i (n) | 


h.(n)  -  h(n)| 


(3.15) 


»*I 

(3.16) 

nel 
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and  therefore 

|hk(n)  -  hk(n)|  >  |hk(n)  -  hk+1(n)|  •  0  nil  (3.17a) 

and 

|hk(n)  -  hk(n)|  -  |hk(n)  -  hk+1(n)|  nel  (3.17b) 

Summing  (3.17a)  and  (3.17b)  over  all  n,  we  obtain 

Ek  -  l  |hk(n)  -  hk(n)|2  >  J  |hfc(n)  -  hk+1(n)|2  (3.18) 


Error  Reduction  in  the  Frequency  Domain 

From  the  triangle  inequality  for  vector  differences: 


|Hk(w)  -  Hk+i (t*>) j  >  | Hk(a>) |  -  !Hk+1(<*>)| 
Therefore,  we  have  from  Parseval's  Theorem  and  (3.19), 


(3.19) 


Eki  l  l\<«)  '  W’l' 


2u 


1 

^27 


|  Hk  (oi )  -  Hk+-|  (o>)  |  dto 


— IT 
IT 


|M(u)  -  ^k+i  (oj)  |  dw 


tak+l 


(3.20) 


Since  Ek  Is  monotone  decreasing  and  Ek>0,  It  has  a  lower  bound  of 

zero.  Consequently,  it  must  converge  to  a  limit  po1nt[30]. 

Although  we  have  shown  that  the  error  Ek  is  monotone  decreasing,  we 

have  not  shown  that  11m  Ek»0  nor  that  11m  ek(w)»eh(w),  i.e.  that  ek(a>) 

k~  k-*» 
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converges  to  eh(u>)  at  each  frequency  w.  In  spite  of  the  possibility  that 

11m  8.  (u)^e.(u)  often.  In  practice,  this  situation  does  not  occur  when  the 
k-*»  *  n 

constraints  of  section  3.1.1  are  Imposed. 

3.2.2  The  OFT  Realization 

Since  our  Iterative  algorithm  will  be  implemented  on  a  digital  com¬ 
puter,  we  can  compute  a  Fourier  transform  at  only  a  finite  number  of  points. 
In  particular,  we  shall  use  the  discrete  Fourier  transform  (DFT).  It  Is 
Important  then  to  investigate  the  Implications  of  this  Implementation. 

We  first  demonstrate  that  under  certain  conditions  uniformly  spaced 
samples  of  the  magnitude  of  the  Fourier  transform  (I.e.,  samples  correspond¬ 
ing  to  the  DFT)  are  sufficient  for  phase  retrieval. 


3. 2. 2.1  Phase  Retrieval  from  Samples  of  the  Magnitude 
Suppose  we  are  given  the  magnitude  of  the  Fourier  transform  of  h(n) 
and  sufficient  a  priori  knowledge  about  the  sequence  so  that  a  unique  phase 
function  can  be  specified.  We  wish  to  show  that  when  our  a  priori  knowledge 
includes  a  finite  duration  constraint  of  M  points,  M  or  more  uniformly 
spaced  samples  of  |H(w)|  at  a»k  *  2irk/N,  ¥wk  in  the  Interval  [0,*],  are  suf¬ 
ficient  to  recover  the  phase  of  the  Fourier  transform  of  h(n). 

To  derive  this  result,  we  examine  the  autocorrelation  function 
R(n)  *  h(n)*h(-n)  whose  Fourier  transform  is  given  by  |H(w)|  .  Since  h(n) 

Is  of  length  M,  the  symmetric  function  R(n)  Is  of  length  2M-1.  At  uni¬ 
formly  spaced  samples  of  the  squared  magnitude,  we  have 


|H(u,k) 


2 


a>k-2irk/N 


M 

£  R(n)  exp£j2irkn/N] 
n«-M 


(3.21) 
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When  N  >_  2M-1,  (3.21)  is  simply  the  DFT  of  R(n).  Therefore,  2M-1  or  more 

2 

samples  of  |H(u)|  In  the  Interval  [0,  2ir]  uniquely  characterizes  both  R(n) 

and  |H(<d)|  .  From  symmetry,  this  condition  is  equivalent  to  specifying  H 

or  more  samples  in  the  Interval  [0,ir][22]. 

Consider  now  two  M  point  sequences  h^(n)  and  h2(n).  Furthermore, 

suppose  that  they  are  recoverable  from  their  respective  magnitude  functions 

and  values  of  the  sequences  over  a  set  I  (e.g.,  the  constraints  of  section 

3.1.1),  and  that  these  values  are  Identical  for  each  sequence:  h^ (n )  - 

h2(n)  *  h(n)  ,  nel.  Suppose  also  that  their  magnitude  functions  are  equal 

at  M  or  more  uniformly  spaced  points  a>k  *  2-irk/N  In  the  interval  [0,ir]. 

Since  both  (n)  and  R2(n)  are  2M-1  points  in  extent  and  are  recoverable 

from  the  given  magnitude  samples.  It  follows  from  above  that  R-|(n)*R2(n) 

2  2 

and  that  |H^(u)|  ■  |H2(w)|  ,  V-^  .  But  we  have  assumed  both  h^n)  and  h2(n) 

are  recoverable  from  their  magnitude  functions  and  the  values  h(n)  for  nel. 
Therefore,  h-j(n)  and  h2(n)  must  be  identical. 

That  is,  contrary  to  our  original  assumption,  there  can  exist  only 
one  M  point  sequence  which  satisfies  the  given  a  priori  knowledge  In  the 
time  domain,  and  has  a  magnitude  function  with  specified  values  at  M  or 
more  uniformly  spaced  points  In  the  interval  [0,ir]. 

3. 2. 2. 2  Iteration  Based  on  Samples  of  the  Magnitude 
From  the  results  of  the  previous  section,  we  are  now  In  a  position  to 
formulate  the  Iterative  procedure  of  section  3.2.1  in  terms  of  a  DFT  reali¬ 
zation.  The  error  function  Is  defined  as  the  mean  squared  difference  be¬ 
tween  samples  of  the  known  magnitude  and  samples  of  the  estimate: 

N-l 

Ek"f  l  |MU)-MkU)|2 

t*o 


(3.22) 
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where  N  Is  the  OFT  length.  h(n)  Is  assumed  causal  and  of  duration  M,  so 
that  N  >  2M-1 . 

Following  steps  (3.14)  through  (3.20)  and  with  the  application  of 
forward  and  inverse  DFTs  we  obtain, 

N-1 

£k  ■  f  lQ  |M(t)expCje(c(£)]  -  Mk(£}expCjek(Jt)]|2 

a  1  V  |Hk(t)  -  HkU)|2 
N  £-0 

N-1 

*  l  |hR(n)  -  hk(n)|2 
n*0 

N-1 

-  JL  *hk(n)  *  Wn)|2 

iHk(£)  -  Hk+1(£)|2 

|M(t)  -  Mk+1(£)|2 

(3.23) 


N-1 

N-1 

-  IT 

l 

£®0 

“Ek+1 

Thus,  as  shown  for  the  continuous  counterpart,  Ek  >_  Ek+j  . 

3.2.3  Examples 

In  this  section  we  Investigate  two  examples  which  Illustrate  the  Iter¬ 
ative  algorithm  for  recovering  phase  from  magnitude.  The  Initial  phase 
estimate  0q(w)  in  both  cases  Is  set  to  zero. 
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Example  3.1 

Consider  a  causal  Infinite  length  mixed  phase  sequence  with  h(0)f0. 

The  first  M+1  points  of  the  sequence  are  constrained  where  M  Is  the  number 
of  zeros  of  H(z)  which  consists  of  two  complex  pole  pairs  at  292  Hz  and 
3500  Hz  and  one  maxlmun  phase  complex  zero  pair  at  2000  Hz.  Consequently, 
the  first  three  points  of  h(n)  are  necessary. 

Since  a  finite  length  sequence  Is  required,  we  assume  h(n)»0  for 
n>256  and  the  OFT  length  Is  set  at  512  points.  This  approximation  does  not 
noticeably  alter  H(u).  Nevertheless,  the  magnitude  of  the  Fourier  transform 
of  the  truncated  sequence  was  Incorporated  in  the  iteration. 

The  sequence  of  functions  log[Mk(u»)]  and  ek(w)  are  depicted  In  Fig. 

3.2  superimposed  on  the  originals  for  5,  10,  20  and  50  Iterations. 
log[MkU)]  and  ek(<i>)  are  Indistinguishable  from  the  originals  after  50 
Iterations. 

Example  3.2 

Consider  a  causal  mixed  phase  finite  length  sequence  where  only  the 
first  point  of  h(n)  Is  constrained  and  the  zero  geometry  of  H(z)  satisfies 
the  condition  of  (3.7).  The  original  sequence  Is  eight  points  In  duration, 
and  a  512  point  OFT  is  used.  The  sequence  of  functions  1  og [Mk (o» ) 3  and  ekU) 
are  depicted  in  Fig.  3.3,  superimposed  on  the  corresponding  originals  for 
1,  2,  4,  10,  80  and  200  Iterations.  log[Mk(w)]  and  0k(u>)  are  Indistinguish¬ 
able  from  the  originals  after  200  iterations. 

3.3  Phase  Estimation  From  Magnitude  by  Transformation 

In  the  previous  sections  we  presented  methods  of  retrieving  the  phase 
from  the  magnitude  function  with  sufficient  a  priori  Information.  An  alter¬ 
native  approach  which  requires  only  a  magnitude  function  Is  to  convert  a 
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(a) 


(b) 


*F1g.  3.2  (a)  Convergence  of  log[Mk(<u)]  in 

example  3.1 ,  (b)  Convergence  of 
ek(u)  in  example  3.1. 


T 


logmagnltude  and  phase  functions  are  measured  in  decibels  and  radians, 
respectively. 


ORIGINAL  (200)4  /  { 


Fig.  3.3  (a)  Convergence  of  log[MkU)]  in 
example  3.2,  (b)  Convergence  of 
Mu>)  in  example  3.2. 
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phase  estimation  problem  to  a  magnitude  estimation  problem  by  transforming 
the  desired  sequence  to  take  on  a  minimum  phase  characteristic.  This 
approach  does  not  require  the  strict  constraints  given  in  section  3.1,  but 
does  require  some  general  knowledge  about  the  sequence  or  its  Fourier 
transform. 

In  the  context  of  deconvolution  of  a  quasi -periodic  waveform,  we  shall 
see  that  our  class  of  transformations  can  be  applied  to  modify  a  desired 
mixed  phase  sequence  while  preserving  the  convolutional  characteristic  of 
the  quasi-periodic  waveform.  After  obtaining  an  estimate  of  the  magnitude 
of  the  modified  (minimum  phase)  sequence,  we  can  obtain  an  estimate  of  its 
phase  by  applying  a  Hilbert  transform.  Performing,  finally,  the  inverse 
operation  to  the  original  transformation  yields  an  estimate  of  the  desired 
sequence,  and  hence  an  estimate  of  its  phase. 

3.3.1  Conditions  for  a  Minimum  Phase  Property 

It  is  desirable  to  detect  a  minimum  phase  characteristic  by  some 
simple  operation  on  the  sequence  or  its  Fourier  or  z-transform  so  that 
appropriate  modifications  for  minimum  phase  conversion  are  easily  derived. 
One  method  of  testing  whether  a  sequence  with  a  rational  z-transform  is 
minimum  phase  is  to  find  the  locations  of  its  poles  and  zeros.  Fortunately, 
a  nunber  of  tests  are  available  that  enable  us  to  determine  whether  or  not 
the  roots  of  the  numerator  and  denominator  polynomials  of  the  z-transform 
lie  within  the  unit  circle.  For  example.  Jury's  criterion[29],  a  counter¬ 
part  of  the  Routh-Hurwitz  criterion[4]  for  stability  of  continuous-time 
systems,  provides  information  on  the  whereabouts  of  the  roots  of  a  poly¬ 
nomial  with  respect  to  the  unit  circle  in  the  z-plane.  Such  tests  enable 
one  to  make  the  denominator  polynomial  of  a  system  function  minimun  phase 
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(i.e.,  stable)  by  appropriately  applying  feedback  to  the  given  system. 

For  our  purposes,  however,  the  discrete-time  version  of  the  Nyquist 
criterion  for  stability  will  be  more  useful  since  It  provides  a  graphical 
interpretation  of  a  necessary  and  sufficient  condition  for  a  minimum  phase 
property.  The  Nyquist  criterion  will  also  be  useful  in  deriving  a  sufficient 
condition  for  a  minimum  phase  z-transform  denoted  as  the  "positive  real" 
constraint.  The  positive  real  constraint  is  then  used  to  derive  other 
sufficient  conditions,  and  also  transformations  to  ensure  the  minimum  phase 
property. 

3. 3. 1.1  The  Nyquist  Criterion 

The  Nyquist  criterion  is  based  on  a  mapping  theorem  by  Cauchy[29].  If 
a  complex  variable  z  in  the  z-plane  describes  a  closed  contour  C1  in  a  pos¬ 
itive  sense,  then  H(z)  will  describe  a  closed  contour  C2  in  the  H(z)  plane. 
The  contour  C2,  the  polar  plot,  will  encircle  the  origin  M  times  in  the 
positive  direction,  where  M  is  the  difference  between  the  number  of  zeros 
and  poles  of  a  rational  z-transform  H(z)  enclosed  by  C-j.  C-j  is  taken  to  be 
the  closed  contour  depicted  in  Fig.  3.4a  where  the  inner  radius  is  the  unit 
circle  and  the  outer  radius  R  is  made  to  approach  infinity.  The  only  con¬ 
tribution  to  the  contour  C2,  i.e.  the  polar  plot,  results  from  the  component 

of  C-,  along  the  unit  circle[29].  The  criterion  requires  that  lim  H(z)  = 

z-*» 

constant  which  holds  In  our  case  since  h(n)  Is  assumed  to  be  causal.  When 
h(n)  Is  also  stable  all  poles  of  H(z)  lie  within  the  unit  circle  and  thus 
the  polar  plot  will  encircle  the  origin  M  times,  where  M  is  the  number  of 
zeros  outside  the  unit  circle. 

A  simple  example  will  be  used  to  Illustrate  the  Nyquist  criterion. 
Consider  a  three  point  sequence  with  z-transform  given  by 


Fig.  3,4  (a)  Contour  within  the  z-plane  required  by 
the  Nyqulst  criterion,  (b)  Polar  plot  of 
H(z)  in  (3.24). 
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H(z)  »  h(0)  +  3.9Z-1  +  4,0z"2  (3.24) 

where  h (0)  is  variable.  For  h(0)*l,  the  z-transform  consists  of  two  zeros 
outside  the  unit  circle  and  two  poles  at  the  origin.  The  polar  plot  of 
H(z)  for  z*exp[jo>]  with  0  <_  u  <_  Zit  is  given  In  Fig.  3.4b.  There  exist  two 
clockwise  encirclements  of  the  origin.  Indicating  the  presence  of  two  zeros 
outside  the  unit  circle. 

One  method  of  modifying  this  sequence  so  that  it  is  transformed  to  a 
minimum  phase  sequence  is  to  increase  h(0)  so  that  there  occur  no  encircle¬ 
ments  of  the  origin.  This  is  equivalent  to  adding  a  positive  constant  to 
H(z)  or  the  real  component  of  H(u>),  Hr(m);  thus  the  entire  polar  plot  is 
shifted  to  the  right  in  the  z-plane.  From  Fig.  3.4b  we  observe  that  the 
leftmost  real  axis  crossing  ( 1 . e.  H- (o»)“0)  occurs  when  Hr(u>)  5  -  3.0.  In 
this  particular  example,  then,  augmenting  h (0 )  so  that  h (0)  >  4  will  ensure 
no  encirclements  of  the  origin.  In  general,  however,  constraining  the 
leftmost  real  axis  crossing  of  the  polar  plot  to  fall  to  the  right  of  the 
origin  is  only  sufficient  for  the  minimum  phase  condition^]. 

Suppose  now  that  h(0)  is  augmented  further  so  that  Hr(u>)  >  0 
Clearly  from  Fig.  3.4b  this  constraint  represents  a  sufficient  condition  on 
Hp (cj)  to  guarantee  no  encirclements  of  the  origin;  that  is  Hr(u>)  must  take 
on  a  negative  value  for  some  w  for  an  encirclement  of  the  origin  to  occur. 

In  our  example,  h (0 )  >  5  guarantees  Hr(w)  >  0  . 

In  summary,  we  can  enumerate  three  Increasingly  restrictive  constraints 
on  H(u>)  for  H(z)  to  be  a  minimum  phase  z-transform  when  h(n)  is  assumed 
causal  and  stable.  We  shall  show  that  these  constraints  imply  increasingly 
restrictive  limitations  as  well  on  the  unwrapped  phase  function  of  H(w). 

The  first  is  the  Nyquist  criterion  which  is  a  necessary  and  sufficient  con- 
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dltion  on  the  polar  plot  of  H(w),  The  polar  plot  can  be  Interpreted  as 
a  tracing  of  the  path  of  H(u)  as  a  function  of  magnitude  and  unwrapped 
phase.  Clearly,  the  unwrapped  phase  function  which  corresponds  to  this 
condition  generally  is  unbounded.  For  example,  the  polar  plot  may  spiral 
outward  and  then  by  necessity  spiral  Inward  without  a  net  ?nc1rclement 
of  the  origin.  The  return  through  an  inward  spiral  must  occur  since  the 
final  value  of  the  unwrapped  phase  equals  its  initial  value  over  the  Inter¬ 
val  [0,ir]  when  the  zeros  of  H(z)  lie  within  the  unit  circle. 

The  second  condition  is  more  restrictive,  follows  directly  from  the 
Nyquist  criterion,  and  is  only  a  sufficient  condition.  This  condition  re¬ 
quires  that  the  leftmost  zero  crossing  of  the  polar  plot  must  fall  to  the 
right  of  the  origin.  As  implied  by  Fig.  3.4b  this  second  condition  guaran¬ 
tees  that  the  unwrapped  phase  of  H(a»)  (assuming  H(0)  >  0  and  thus  0h(O)=*O) 
cannot  exceed  ir  in  absolute  value:  |eh(ui)|  <  it  . 

The  third  condition  is  somewhat  more  restrictive  than  the  second, 
again  follows  directly  from  the  Nyquist  criterion,  and  requires  that  the 
real  component  of  H(u)  be  positive.  As  implied  by  Fig.  3.4b  this  positivity 
constraint  guarantees  that  the  unwrapped  phase  cannot  exceed  ir/2 :|b^(c*») |<tt/2. 
We  shall  see  in  the  next  section  that  the  implications  of  this  positivity 
constraint  are  even  greater  since  H(z)  under  this  condition  falls  into  a 
class  denoted  as  "positive  real  z- transforms". 

The  phase  restrictions  due  to  the  latter  two  conditions  guarantee  that 
the  unwrapped  phase  equals  the  principal  value  of  the  phase.  Consequently, 
the  unwrapped  phase  may  be  computed  at  samples  by  simply  an  arctangent 
routine. 
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3. 3.1. 2  The  Positive  Real  Constraint 

A  rational  function  of  a  complex  variable  s,  F(s)  which  is  real  for 
real  values  of  s,  and  whose  real  part  is  positive  for  all  values  of  s  with 
a  positive  real  part,  is  called  a  positive  real  function.  Functions  of 
this  sort  play  an  important  part  in  electrical  network  theory  and  have 
been  studied  extens1vely[9].  In  this  section  we  first  briefly  review  the 
properties  of  a  positive  real  function.  The  minimum  phase  characteristic 
of  such  a  function  is  proven  in  a  new  way  through  the  continuous-time  ver¬ 
sion  of  the  Nyquist  criterion.  We  then  proceed  to  develop  an  analogous 
theory  for  complex  functions  with  respect  to  the  unit  circle.  Specifically, 
a  positive  real  rational  z-transform  for  a  causal  and  stable  sequence  Is 
defined  by  the  following  two  properties: 

(i)  H(z)  is  real  for  real  values  of  z 
(1  i )  Hr(z)  >  0  for  |z|  _>  1 

The  first  property  Implies  simply  that  the  sequence  is  real.  We  saw  in  the 
previous  section  that  when  h(n)  is  causal  and  stable  and  Hr(u))  >  0,  h(n) 
is  a  minimum  phase  sequence.  We  shall  show  that  these  conditions  are  nec¬ 
essary  and  sufficient  for  H(z)  to  be  a  positive  real  z-transform.  Thus  in¬ 
vestigation  of  the  positivity  of  Hr(z)  over  the  entire  region  |z|  >.  1  is 
not  necessary  in  determining  whether  H(z)  is  positive  real. 

Consider  now  the  positive  real  function  F(s)  whose  properties  were 
described  above.  Gulllemln  has  shown  that  a  positive  real  function  cannot 
have  poles  in  the  right  half  s-plane  (RHP)  since  the  presence  of  a  pole 
Implies  fj,(s)  <  0  for  some  s  in  the  region  of  the  pole[9].  Clearly  F(s) 
cannot  have  zeros  in  the  RHP  and  therefore  F(s)  represents  a  minimun  phase 
function.  Through  appropriate  use  of  the  theory  of  functions  of  a  complex 
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variable.  It  Is  possible  to  show  that  a  necessary  and  sufficient  condition 
for  F(s)  to  be  positive  real  is  that  F(s)  is  positive  on  the  ju  axis  and 
analytic  for  Re(s)  >  0[9],  Therefore,  an  alternative  way  of  proving  the 
minimum  phase  characteristic  of  positive  real  functions  is  to  apply  the 
continuous- time  version  of  the  Nyquist  criterion  to  the  polar  plot  of  F(s) 
along  the  ju>  axis.  The  analytlcity  constraint  for  Re(s)  >  0  guarantees 
no  poles  in  the  RHP.  The  positivity  constraint  guarantees  that  there  can 
exist  no  encirclements  of  the  origin  by  the  polar  plot,  C2, 
where  C-j  is  taken  to  be  the  boundary  of  the  RHP,  Including  the  ju  axis. 
Therefore,  the  RHP  Is  also  zero  free  and  F(s)  Is  a  minimum  phase  function. 

The  discrete-time  counterpart  to  this  theory  easily  follows: 


A  z-transform  H(z)  Is  positive  real  with  respect  to 
the  unit  circle  if  and  only  If  H(z)  is  analytic  for 
|z|  >_  1  and  Hp(z)  >  0  on  the  unit  circle. 


The  proof  of  the  "only  if"  component  of  this  proposition  is  Identical 
to  that  byGuillemin  for  a  positive  real  function[9],  and  will  not  be  pre¬ 
sented  because  of  the  little  insight  gained.  The  proof  of  the  "If"  com¬ 
ponent,  on  the  other  hand,  provides  a  rather  Interesting  view  Into  the  nat¬ 
ure  of  a  positive  real  z-transform.  The  proof  Is  also  similar  to  that  by 
Gulllemln,  but  is  modified  to  address  the  discrete-time  nature  of  the  prob¬ 
lem. 

The  proof  calls  upon  the  Extremum  Theorem[9]  of  complex  analysis 
which  concerns  the  real  part  of  a  function,  analytic  over  a  given  region. 
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and  on  the  boundary  of  that  region.  The  theorem  states  that  the  largest 
and  smallest  values  which  the  real  part  (or  Imaginary  part)  assumes  through¬ 
out  the  given  region  Including  Its  boundary  must  lie  on  the  boundary. 
Therefore  the  smallest  values  which  the  real  part  assumes  throughout  this 
region  Including  Its  boundary  must  lie  on  the  boundary.  In  our  particular 
problem  the  smallest  value  assured  by  the  real  part  of  H(z)  In  the  Infin¬ 
itely  large  annulus  whose  boundary  includes  the  unit  circle  and  a  circle 
at  Infinity  must  occur  then  on  this  boundary.  Since  Hr(z)  >  0  on  the  unit 
circle,  the  value  of  H(z)  at  Infinity  Is  also  positive: 

IT 

lim  H(z)  -  h (0 )  *  H  (u>)  da>  >  0  (3.25) 

-IT 

Thus  our  entire  boundary  Is  positive  and  from  the  Extremum  Theorem  no  point 
within  this  boundary  can  be  zero  or  negative. 

Clearly,  we  can  use  the  positivity  of  Hp(u>)  and  analyticity  of  H(z) 
for  | z [  _>  1  to  prove  that  a  positive  real  z-transform  is  minimum  phase.  We 
have  in  fact  already  shown  this  property  in  the  preceding  section  through 
the  Nyquist  criterion. 

3.3. 1.3  A  Sufficient  Condition  on  h (0) 

The  positive  real  constraint  is  useful  in  testing  for  a  minimum  phase 
characteristic.  It  is  also  useful  In  deriving  other  sufficient  conditions 
for  ensuring  this  property.  In  particular,  we  shall  capitalize  on  the 
positive  real  constraint  In  proving  the  following  sufficient  condition: 

A  causal  sequence  h(n)  Is  minimum  phase  if  the  absolute  value  of  h(0) 

Is  greater  than  the  sum  of  the  absolute  value  of  the  remaining  terms: 

|h(0)|  >  I  |h(n) | 
nM 


(3.26) 
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We  shall  show  by  example  that  this  condition  Is  more  restrictive  than  the 
positive  real  constraint. 

Since  h(n)  is  causal  and  stable,  H(z)  is  analytic  for  |z|  _>  1 .  Without 
loss  of  generality  we  also  assume  that  h(0)  is  positive.  Suppose  now  that 
for  some  &>,  the  positive  real  condition  doesn't  hold: 


Therefore , 


H„(<i))  ■  h (0 )  +  T  h(n)  cosun  <  0 
r  n>l 


T  h(n)  cos un  <  -  h(0) 
n>l 

Since  h (0 )  >  0,  we  have 


(3.27) 


(3.28) 


I  y  h(n)  cosun  |  >  |h(0)|  (3.29) 

n>l 

and  so,  since  [cosunj  <_  1 


l  !h(n)|  i  l  |h(")  cosu>n| 
n^l  nH 


>_  l  h(n)  cosun 
n>_l 


I  h (0 )  | 


(3.30) 


In  conclusion, 

| h (0) |  <  l  |h(n) |  (3.31) 

n>l 

which  is  a  contradiction  to  our  original  assumption.  Therefore,  Hr(a»)  >  0, 
^  and  since  H(z)  Is  analytic  for  |z|  >  1.  H(z)  Is  positive  real,  and  con- 

U) 

sequently  minimum  phase. 

The  positive  real  constraint,  however,  in  general  does  not  imply  con- 
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dltlon  (3.26).  A  counterexample  Is  the  sequence  h(n)  *  an  u(n)  where 
.5  <  a  <  1  and  where  u(n)  Is  the  unit-step  sequence.  Summing,  we  obtain 

l  |h(n)|  *  l  a" 
n>l  n*l 

■  a/d-a) 

>  |h(0)|  (3.32) 

for  .5  <  a  <  1.  But,  Hr(w)  >0f 

3.3.2  Transformations  for  a  Minimum  Phase  Sequence 
In  this  section  we  present  a  number  of  techniques  for  transforming  an 
arbitrary,  causal  mixed  phase  sequence  to  a  minimum  phase  sequence.  These 
transformations  are  based  on  the  conditions  discussed  In  the  previous  sec¬ 
tion,  and  are  restricted  to  a  class  which  has  two  specific  properties: 

(1)  A  transformation  T  must  be  invertible 

(ii)  In  the  context  of  deconvolution  of  a  quasi -periodic  waveform, 

T  must  modify  the  desired  sequence,  while  preserving  the  con¬ 
volutional  characteristic  of  the  original  waveform.  That  is, 
a  transformation  T  must  be  such  that 

x(n)  »  T[x(n)] 

■  T[h(n)*p(n)] 

-  T[h(n)]*p(n)  (3.33) 

where  T'1  exists,  l.e. ,  T_1[T[h(n)]]  »  h(n),  and  p(n)  may  or  may  not  equal 

/\ 

p(n),  but  consists  of  a  train  of  equally  spaced  samples  with  spacing  P. 
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We  shall  consider  two  such  transformations: 

(i)  Exponential  weighting 

(11)  Addition  of  a  "reference"  signal 

3. 3.2.1  Exponential  Weighting 

Consider  weighting  a  mixed  phase  sequence  with  a  decaying  exponential 

[32]: 

h(n)  *  anh(n)  (3.34) 

This  operation  results  in  preservation  of  the  convolutional  property  and  in 
individual  weighting  of  both  components  of  a  quasi -periodic  waveform: 

x(n)  *  on(h(n)*p(n)) 

=  txnh(n)*anp(n)  (3.35) 

The  z-transform  of  h(n)  *  anh(n)  is  given  by 

H(z)  =  l  anh(n)z“n  »  H(a_1z)  (3.36) 

n 

Therefore  if  H(z)  has  a  pole  or  zero  at  z*zQ,  H(z)  has  a  pole  or  zero  at 
azQ.  h(n)  can  then  be  made  minimum  phase  by  multiplying  x(n)  by  a11  where 
a  is  small  enough  to  move  the  pole  or  zero  of  h(n)  with  the  greatest  magni¬ 
tude  inside  the  unit  circle:  (n)*onh(n).  h(n)  can  always  be  recovered 

from  hmp(n)  by  inverse  transformation:  h(n)“a“nh[np(n). 

One  drawback  to  this  technique  is  that  often  the  required  a  is  small 
enough  to  cause  trouble  with  rounding  error  when  implemented  with  a  digital 
computer[32],  A  second  drawback  is  that  o  depends  on  knowledge  of  zero 
locations. 
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3. 3. 2. 2  Addition  of  a  Reference  Signal 

A  different  approach  is  to  add  a  causal  "reference"  signal  g(n)  to 
h(n)  to  create  a  minimum  phase  sum  h(n)  =  g(n)  +  h(n).  Such  transforma¬ 
tions  are  also  suitable  for  preserving  the  features  of  a  quasi -periodic 
waveform  since  with  a  priori  knowledge  of  p(n),  we  can  form  the  sequence 

x(n)  *  h(n)*p(n)  +  g(n)*p(n) 

a  (h(n)  +  g(n))*p(n)  (3.37) 

One  simple  way  of  guaranteeing  that  h(n)  is  minimum  phase  is  to  ensure 
that  the  polar  plot  of  HU)  has  no  encirclements  of  the  origin  by  imposing 
the  constraint  HrU)  >0.  An  alternative  is  to  force  the  leftmost  real  axis 
crossing  of  the  polar  plot  to  lie  to  the  right  of  the  origin  in  the  z-plane. 
The  latter  condition  can  be  checked  by  Investigation  of  the  gain/phase 
margin  for  stability  and  thus  directly  from  a  logmagnltude  and  unwrapped 
phase  function [4], 

A  third  approach  is  to  add  a  minimum  phase  signal.  We  shall  show  with 
the  aid  of  the  positive  real  condition  that  when  G(z)  is  minimum  phase  and 
|GU)|  >  | H (oj )  |  ,  ¥■  ,  H(z)  *  G(z)  +  H(z)  is  minimum  phase.  To  see  this,  we 
express  H(z)  by 

H(z)  -  G(z)(l  +  H(z)/G(z) )  (3.38) 

G(z)  Is  minimum  phase,  so  the  remaining  problem  is  to  show  that  (1+H(z)/G(z)) 
Is  minimum  phase.  Since  |HU)/G(u>)|  <  1 ,  the  polar  plot  of  HU)/GU)  must 
lie  within  the  unit  circle  |z|  <  1.  Therefore,  l+HU)/G(u>)  lies  to  the 
right  of  the  imaginary  axis  In  the  z-plane,  so  that  Re(l+HU)/GU))  >  0. 

Now,  both  H(z)  and  G(z)  represent  causal  sequences  and  since  G(z)  is  mini¬ 
mum  phase,  G"^(z)  also  represents  a  causal  sequence  (see  section  2.2.1). 
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Therefore  l+H(z)/G(Z)  corresponds  to  a  causal  sequence,  and  must  be  analytic 
for  \z\  _>  1.  Consequently,  1+H(z)/G(z)  Is  positive  real,  and  thus  minimum 
phase. 

We  have  already  seen  an  example  in  section  3. 3. 1.1  of  transforming  a 
mixed  phase  sequence  to  a  minimum  phase  sequence  by  augmenting  h(0)  with  a 
positive  constant:  g(n)  »  A6(n).  The  conditions  we  have  been  discussing 
together  with  that  of  section  3. 3. 1.3  represent  restrictions  on  the  values 
of  A  of  different  degrees. 

As  a  point  of  interest,  we  can  show  that  modifying  the  value  of  h(0) 
to  form  a  minimum  phase  sequence  is  analogous  to  changing  the  gain  of  a 
feedback  control  system  for  stability.  Factoring  out  A  in  h(n)*A6(n)+h(n) , 
we  obtain 

H(z)  ■  A(1  +  A~ 1 H ( Z ) )  (3.39) 

A"^H(z)  is  reminiscent  of  an  "open  loop  transfer  function"  where  A"^  is 


the  "open  loop  gain"[4]. 
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CHAPTER  4 

MIXED  PHASE  DECONVOLUTION 

The  problem  of  deconvolving  a  sequence  from  a  quasi -periodic  waveform 
was  described  in  section  2.1.  In  the  frequency  domain  the  problem  of  decon¬ 
volution  is  twofold:  to  estimate  the  magnitude  and  phase  of  a  Fourier 
transform.  In  Chapter  3  we  assumed  that  the  desired  magnitude  is  known  or 
measurable  and  that,  with  a  priori  information  about  the  sequence  or  phase 
function,  we  found  that  unambiguous  phase  recovery  is  possible.  We  also 
investigated  methods  of  creating  a  situation  where  the  phase  can  be  recovered 
from  the  magnitude.  That  is,  we  derived  transformations  which  convert  a 
mixed  phase  sequence  to  a  minimum  phase  sequence,  and  which  are  applicable 
to  magnitude- only  deconvolution. 

One  technique  which  is  potentially  capable  of  directly  and  accurately 
estimating  both  the  magnitude  and  phase  of  a  Fourier  transform  from  a  quasi- 
periodic  waveform  is  homomorphic  deconvolution [24, 32].  A  drawback  to  this 
direct  approach,  however,  is  the  requirement  of  an  unwrapped  phase  function. 
The  unwrapped  phase  is  generally  difficult  to  compute  at  samples  due  to 
modulo  2-n-  considerations.  Available  unwrapping  algorithms  are  prone  to 
error  when  frequency  sampling  is  not  "sufficiently"  dense[32,34]  or  the 
Fourier  transform  contains  regions  of  low  energy  which  are  particularly 
susceptible  to  degradation  by  quantization  noise.  Moreover,  as  we  shall 
see,  the  envelope  (i.e.,  slowly  varying  component)  of  the  unwrapped  phase 
is  quite  sensitive  to  small  changes  of  a  sequence  in  the  time-domain.  This 
sensitivity  leads  to  an  inherently  ill-conditioned  problem. 

As  a  result,  we  wish  to  either  avoid  the  use  of  an  unwrapped  phase 
function,  or  compute  this  function  from  only  frequency  regions  with  a  high 
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signal-to-noise  ratio.  The  latter  approach  was  first  taken  by  Tribolet[33] 
in  confronting  the  problem  of  low-pass  and  high-pass  filtered  seismic  data. 
His  approach  entailed  shifting  and  stretching  the  signal's  passband  to 
occupy  the  entire  frequency  interval  [Q.tt],  while  preserving  the  convolu¬ 
tional  characteristic  of  the  original  sequence.  Our  problem,  however,  is 
not  suitable  to  this  approach  since  the  periodic-like  nature  of  our  wave¬ 
forms  corresponds  in  the  frequency  domain,  to  multiple  bands  of  low  signal- 
to-noise.  In  Chapter  5  we  develop  a  technique  within  this  same  philosophy, 
which  addresses  the  harmonic  structure  of  our  spectra. 

The  alternative  approach  is  to  avoid  the  issue  of  phase  unwrapping  by 
conversion  to  a  magnitude-only  deconvolution  problem.  Applying  the  trans¬ 
formations  of  section  3. 3. 2, 2,  we  can  perform  homomorphic  deconvolution 
without  an  unwrapped  phase.  Note  that  although  we  shall  illustrate  magni¬ 
tude-only  deconvolution  through  homomorphic  filtering,  any  deconvolution 
technique  yielding  a  minimun  or  zero  phase  solution  is  applicable. 

When  such  an  approach  is  difficult  to  take,  as  for  example  with  in¬ 
adequate  knowledge  of  p(n),  we  must  compute  DFT  samples  of  an  unwrapped 
phase  in  directly  estimating  phase  by  homomorphic  deconvolution.  The  need 
of  an  unwrapped  phase  arises  also  in  the  following  chapter  where  a  phase 
estimate  is  deduced  from  only  harmonic  samples  of  a  Fourier  transform. 

The  second  major  area  of  this  chapter  presents  a  new  technique  for 
determining  the  unwrapped  phase  which  does  not  require  conventional  modulo 
2i  considerations.  In  developing  this  technique,  we  first  investigate  mag¬ 
nitude  retrieval  from  the  phase  of  a  Fourier  transform.  In  a  manner  simi¬ 
lar  to  our  development  in  sections  3,1  and  3.2  we  formulate  a  linear  iter¬ 
ative  procedure  which  retrieves  the  magnitude  from  the  phase,  and  which 


» 


L 


constrains  the  values  of  a  sequence  over  a  specified  region. 


Under  a  causality  constraint,  the  algorithm  provides  an  alternative  to 
the  Hilbert  transform  in  obtaining  the  magnitude  from  the  phase  for  a  mini¬ 
mum  phase  sequence.  However,  It  has  the  advantage  of  not  requiring  an  un¬ 
wrapped  phase  function,  but  only  the  principal  value  of  the  phase.  Under  a 
finite  length  constraint,  the  algorithm  is  capable  of  unambiguously  recov¬ 
ering  a  magnitude  from  an  arbitrary  phase. 

Finally,  in  the  context  of  phase  unwrapping,  the  algorithm  serves  as  a 
major  component  within  our  new  phase  unwrapping  algorithm. 

4.1  Phase  Estimation  by  Homomorphic  Deconvolution 

When  we  confront  the  problem  of  filtering  signals  which  have  been  added, 
we  often  use  a  linear  filter.  Extracting  the  phase  or  the  magnitude  of  H(u>), 
on  the  other  hand,  is  a  nonlinear  problem  since  HU)  and  PU)  are  multlplica- 
tively  combined.  The  approach  of  homomorphic  system  analysis  proposed  by 
0ppenheim[24]  transforms  this  nonlinear  problem  to  a  linear  filtering  prob¬ 
lem.  In  this  section  we  first  review  the  theory  of  multiplicative  homomor¬ 
phic  systems,  and  "direct"  phase  estimation  by  homomorphic  filtering.  We 
then  proceed  to  illustrate  the  sensitivity  inherent  in  this  direct  approach 
due  to  the  requirement  of  an  unwrapped  phase  function.  Finally,  a  specific 
transformation  of  section  3. 3. 2. 2  is  applied  in  demonstrating  the  indirect 
approach  of  phase  estimation  by  magnitude-only  homomorphic  deconvolution. 

4.1.1  Theory 

Much  of  the  groundwork  for  multiplicative  homomorphic  systems  has  been 
laid  in  Chapter  2.  The  purpose  of  such  systems  Is  to  transform  the  problem 
of  separating  the  multiplicative  components  of  xU)  3  HU)P(u>)  to  a  problem 
of  linearly  filtering  additive  components.  In  particular,  the  complex  log- 


66 


arithm  provides  the  means  of  obtaining  this  additivity: 

1og[X(z)]  *  log[H(z)]  +  log[P(z)]  (4.1) 

■  log |H(z) |  +  log|P(z) |  +  j(eh(z)  +  ep(z) ) 

The  additivity  of  the  imaginary  components  (i.e.,  the  phase  components) 

of  the  complex  logarithm  of  X(z)  holds  when  the  unwrapped  phase  Is  used  in 

n 

(4.1).  A  linear  phase  component  z  0  may  also  be  included  when  its  unwrapped 
phase  is  defined  to  be  odd,  but  discontinuous  at  ir  .  This  term,  neverthe¬ 
less,  should  be  removed  since  it  interferes  with  the  estimation  of  eh(u>) 
[32]. 

Under  these  conditions,  the  complex  logarithmic  operation  represents  a 
homomorphic  system  which  maps  multiplication  to  addition.  Taking  the  in¬ 
verse  transform  of  (4.1)  we  obtain  the  real  sequence 

x(n)  -  h(n)  +  p(n)  (4.2) 

x(n)  is  termed  the  "complex  cepstrum"  to  emphasize  the  use  of  the  complex 
logarithm.  We  use  the  term  "real  cepstrum"  when  the  real  logarithm  is  used 
and  thus  only  the  magnitude  of  the  Fourier  transform  of  x(n)  is  retained. 
Likewise,  "phase  cepstrum"  refers  to  the  case  where  only  the  phase  is  re¬ 
tained. 

For  a  rational  H(z),  (4.2)  can  be  expressed  in  terms  of  the  complex 
cepstrum  of  the  minimum  and  maximum  phase  components  of  the  desired  se¬ 
quence  and  the  pulse  train: 

x(n)  *  (log[A])i(n)  +  hmin(n)  +  hmax(n)  +  p(n)  (4.3) 

where  hm...(n)  and  h _ (n)  are  normalized  so  that  h,„(0)  *  hi  (0)  »  1. 

min  max  min  max 
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Using  the  well-known  power  series  expansions: 
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and 
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From  (2.19)  and  (2.20),  It  follows  that 
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(4.5a) 


(4.5b) 


where  u(n)  Is  the  unit-step  function  and  |ak|,  |bk|,  |ck|,  and  |dk|  are  less 
than  unity. 

The  properties  of  the  complex  cepstrum  that  are  of  Importance  In  the 
sequel  can  be  sunmarized  as  follows: 


PI  h(n)  decays  at  least  as  fast  as  1/n.  Specifically, 

,n> 


I h(n) |  <  C 


0 


n 


(4.6) 


where  C  is  a  constant  and  e  equals  the  maximum  of  |ak|,  | | , 
|ck|,  and  jdk|. 


P2  The  complex  cepstrum  of  a  minimum  phase  sequence  Is  stable 
and  causal,  and  likewise  the  complex  cepstrun  of  a  maximum 
phase  sequence  is  stable  and  anti-causal. 


P3  The  complex  cepstrum  p(n)  of  a  train  of  equally  spaced 


Suppose  that  h(n)  decays  rather  quickly  and  that  the  sample  spacing 
of  p(n)  is  sufficiently  large  so  that  h(n)  and  p(n)  do  not  considerably 

A  A 

overlap.  Under  these  conditions  we  can  filter  h(n)  from  p(n)  by  applying 
a  "7ow-t1me  gate"  to  x(n).  An  estimate  of  h(n)  is  computed  by  transform¬ 
ing  the  estimate  of  h{n),  exponentiating,  and  inverse  transforming. 
Estimation  of  h(n)  is  equivalent  to  obtaining  an  estimate  of  the 

A 

even  and  odd  components  of  h(n).  These  components  correspond  to  log|H(u>)| 
and  respectively.  Therefore,  we  might  estimate  the  two  components 

independently  with  low-time  gates  tailored  to  each  separately. 

4.1.3  The  Heuristics  of  Unwrapped  Phase  Sensitivity 
Suppose  that  the  sequence  x(n)  *  h(n)*p(n)  is  modified  by  an  additive 
disturbance  d(n): 

x(n)  =  h(n)*p(n)  +  d(n)  (4.7) 

In  the  frequency  domain,  (4.7)  is  given  by 
X(ui)  *  H(aj)P(u)  +  D(u) 

=  H(«)P(a>)(l  +  D(u)/H(«)P(w)) 

■  H(ut)P(w)E(u») 

where , 

EU)  -  (1  +  D(uj)/H(u>)P(u>) ) 

Therefore , 


(4.8) 

(4.9) 


logCx(u)]  *  log[H(u)P(u)]  +  log[E(a))] 


(4.10) 
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Thus,  the  effect  of  D(u)  on  the  complex  cepstrum  is  strongly  dependent 
on  the  ratio  D(a>)/H(u>) P(o») . 

For  homomorphic  deconvolution  to  be  robust  in  the  presence  of  a  dis¬ 
turbance  d(n),  the  slowly  varying  components  (i.e.,  the  envelope)  of  the 
complex  logarithm  should  not  be  susceptible  to  a  large  change  with  a  small 
change  in  x(n)  due  to  d(n).  We  wish  to  preserve  the  slowly  varying  com¬ 
ponents  of  log { X (o>)  |  and  ex(ui)  since  they  map  to  the  low-time  region  of  the 
complex  cepstrum.  In  the  context  of  quasi -periodic  waveforms  the  follow¬ 
ing  observations  were  made: 

(I)  "Small M  disturbances  in  the  sequence  x(n)  tend  to  propagate 
small  changes  In  the  envelope  of  log  |X(u>)|  and  thus  small 
changes  In  the  low-time  region  of  the  real  cepstrum. 

(II)  Small  disturbances  in  x(n)  often  propagate  "large"  changes 
In  the  envelope  of  the  unwrapped  phase  of  0x(<j)  and  thus 
large  changes  in  the  low-time  region  of  the  phase  cepstrum. 

(iii)  Large  changes  in  the  unwrapped  phase  envelope  of  ex(“)  may  or 
may  not  be  mapped  through  low- time  cepstral  gating  to  large 
changes  in  the  impulse  response  estimate[28,35]. 

In  suirmary,  an  estimate  of  h (n )  obtained  by  homomorphic  deconvolution 
which  requires  a  magnitude  function  only  is  less  susceptible  to  degrada¬ 
tion  by  time-domain  disturbances  than  an  estimate  based  on  magnitude  and 
phase. 

The  sensitivity  of  the  unwrapped  phase  function  is  understood  by  the 
following  heuristic  argument.  Consider  a  region  where  |D(w)/H(u))PU) |  <  1. 
Under  this  constraint  the  unwrapped  phase  of  E(w)  cannot  exceed  */2  in 
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absolute  value,  as  depicted  in  Fig.  4.1a,  where  the  "noise  vector", 
D(w)/H(<d)P(u>)  remains  smaller  In  magnitude  than  the  "signal  vector",  1+jO. 

When  the  noise  vector  exceeds  the  signal  vector  in  magnitude,  as  for 
example  In  a  low-energy  region  of  H(w)P(w),  the  polar  plot  of  E(a>)  may 
encircle  the  origin  resulting  in  a  2ir  jump  in  the  unwrapped  phase  of  E(a>). 

Such  2ir  jumps  can  accumulate  and  significantly  distort  the  envelope  of 
the  unwrapped  phase  of  ex(u»)»  and  thus  the  low-time  region  of  the  phase 
cepstrum.  Gating  the  phase  cepstrum  may  therefore  not  extract  an  accurate 
estimate  of  the  desired  phase,  0h(w). 

4.1.4  Magnitude-Only  Deconvolution  by  Transformation: 

_ An  Indirect  Approach _ _ _ 

Consider  modifying  a  causal  sequence  h(n)  by  adding  to  it  a  sequence 
g(n)  so  that  the  resulting  sum  v(n)  ■  h(n)  +  g(n)  is  minimum  phase.  In 
particular,  we  choose  g(n)  *  A6(n).  From  (3.36)  to  preserve  the  convolu¬ 
tional  property  of  x(n),  and  to  appropriately  modify  h(n),  we  add  to  x(n) 
the  sequence  A<s(n)*p(n)  to  obtain: 

y(n)  =  (h (n )  +  A<5(n)}  *  p(n) 

*  v(n)*p(n)  (4.11) 

Since  As(n)  is  such  that  v(n)  is  minimum  phase,  from  property  P2  of 

A 

section  4.1.1,  the  complex  cepstrum  of  v(n),  v(n)  is  causal  and  thus  com¬ 
pletely  characterized  by  its  even  component  or  equivalently  by  log|V(ui)|. 
From  (2.11a)  the  complex  cepstrum  of  v(n)  is  therefore  expressed  by 

v(n)  *  L(n)vr(n)  (4.12) 

where  vf(n)  is  the  real  cepstrum  of  v(n)  and  L(n)  is  the  causal  sequence 
given  by 


Fig,  4,1  (a)  Polar  plot  of  E(u)  in  (4.9)  with  a 

small  disturbance  (no  origin  encirclement) 
(b)  Same  as  (a)  with  a  large  disturbance 
(origin  encirclement) 
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p  n-0 

l(n)  »  i  2  n>0  (4.13) 

ijO  n<0 

Since  H(z)  is  assumed  rational,  V(z)  is  also  rational  and  therefore  v(n) 
satisfies  properly  PI  of  section  4,1.1. 

The  real  cepstrum  of  p(n)  is  the  even  component  of  the  complex  cep- 

A 

strum  of  p(n),  p(n)  and  from  property  P3  of  section  4.1.1  consists  of 
equally  spaced  samples  with  spacing  P.  Consequently,  when  v(n)  decays 
rather  quickly,  the  real  cepstra  vr(n)  and  pp(n)  do  not  significantly  over¬ 
lap.  From  this  property  and  (4.12)  we  can  obtain  an  estimate  of  v(n)  by 
gating  the  real  cepstrum  of  y(n): 

A  A 

v(n)  =  g(n)  yr(n) 

A  A 

*  g(n)vp(n)  +  g(n)pp(n) 
f  v(n)  0  <  n  <  P 

=  1  v(0)  +  p(0)  n*0  (4.14) 

v.  0  n<0,  n>p 

where 

fl  n=0 

g(n)»<2  0  <  n  <  P  (4.15) 

(.0  n<0,  n>P 

A 

With  p(0)*0,  we  then  obtain  an  estimate  of  v(n)  by  transforming 
9(n)yr(n),  exponentiating  and  inverse  transforming.  Finally,  an  estimate 
of  h(n)  is  obtained  through  the  inverse  operation  of  subtracting  A<s(n) 
from  the  estimate  of  v(n). 
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With  respect  to  the  choice  of  the  value  A,  a  number  of  conments  are 
In  order.  We  saw  In  section  3. 3.2.2  that  the  value  of  A  may  be  chosen  to 
satisfy  any  one  of  a  number  of  sufficient  conditions  for  ensuring  a  mini¬ 
mum  phase  characteristic.  The  extreme  case  of  choosing  A  so  that  the 
Nyqulst  criterion  Is  just  met  should  be  avoided.  This  Is  because  such  a 
choice  places  a  zero  of  V(z)  "just  on"  the  unit  circle.  From  property  PI 
of  section  4.1.1,  the  complex  cepstrum  of  v(n)  decays  as  approximately 
C|p-|  ,  and  so  the  requirement  that  v(n)  decay  rather  quickly  Is  not  satis 
fled. 

The  alternative  extreme  is  to  allow  the  value  of  A  to  become  very 
large.  With  this  choice,  the  complex  cepstrum  of  y(n)  has  an  interesting 
and  rather  useful  property.  Replacing  az"1  by  H(z)/A  in  the  logarithmic 
expansion  of  (4.4a),  and  assuming  |H(z)/A|  «  1,  we  obtain 

log[Y(z)]  »  log [A  +  H(z)]  +  log[P(z)] 

*  log[A]  +  log[l  +  H(z)/A]  +  log[P(z)] 


’  log[A]  +  H(z)/A  +  log[P(z)] 


(4.16) 


y(n)  is  then  written  approximately  as 

’  log[A]  +  h(0)/A  +  p(0) 
y(n)  =  <  h(n)/A  +  p(n) 

.  p(n) 


n»0 

n>0 

n<0 


(4.17) 


Therefore,  h(n)  can  be  recovered  "almost  exactly"  and  directly  from  the 
complex  cepstrum  by  setting  y(n)*0  for  n»kP,  where  k»+l,+2,  ...,  subtract 
ing  log[A]  +  p(0)  from  y(0),  and  scaling  the  result  by  A.  An  equivalent 
operation  can  also  be  carried  out  on  the  real  cepstrum. 
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A 

In  order  to  eliminate  the  need  of  a  priori  knowledge  of  p(0)  in  the 
above  algorithms,  we  apply  g(n)  a  A<s(n)  to  a  version  of  h(n)  which  is 
shifted  to  the  right: 

v(n)  =  A<s(n)  +  h(n-nQ)  (4.18) 

For  this  case,  we  can  show  that  the  estimate  obtained  from  (4.14) 

A 

for  n>nQ  is  given  approximately  by  exp[p(0)]  h(n-nQ).  Thus,  only  a  scaling 
degradation  occurs  so  that  the  desired  phase  is  preserved. 

4.1.5  Examples 

We  now  consider  two  examples  of  the  techniques  discussed  within  this 
section.  We  shall  compare  estimates  from  the  direct  and  indirect  approaches 
of  sections  (4.1.2)  and  (4.1.4)  obtained  by  gating  the  complex  and  real 
cepstrum,  respectively,  with  and  without,  a  white  noise  disturbance. 

Example  4.1 

Consider  H(z)  with  two  complex  pole  pairs  at  292  Hz  and  3500  Hz,  and 
a  maximum  phase  complex  zero  pair  at  2000  Hz.  p(n),  the  pulse  train  is 
given  by 

p(n)  »  ,3<5(n)  +  5(n-64)  +  ,55(n-128)  (4.19) 

where  P=64.  In  order  to  reduce  the  effect  of  si  delobe  interference  of 
the  Fourier  transform  of  a  rectangular  gate,  the  complex 
cepstrum  Is  multiplied  by  a  Hamming  window  expressed  by 

f.54  -  .46  cos  C 60  £  "  1  60 
w(n)  =  < 


0 


|n|  >  60 


(4.20) 
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The  original  continuous  unwrapped  phase  and  its  estimate  derived  from 
the  direct  approach  are  Illustrated  in  Figs.  4.2a  and  4.2b,  respectively. 

Next,  the  indirect  approach  Is  taken  where  p(n)  is  appropriately 

A 

scaled  so  that  p(0)a0.  We  set  A*. 05  to  create  a  minimum  phase  sequence, 
v(n)aAs(n)  +  h(n).  The  real  cepstrum  of  x(n)  is  multiplied  by  w'(n)*w(n)L (n) 
where  L(n)  is  given  in  (4.13)  and  w(n)  is  given  in  (4.20).  The  resulting 
continuous  phase  estimate  is  depicted  in  Fig.  4.2c  and  is  almost  indistin¬ 
guishable  from  that  derived  from  the  direct  approach. 

Example  4.2 

This  example  Is  identical  to  example  4.1,  but  now  we  add  white  noise 
to  x(n)  with  a  40  db  S/N.  The  original  unwrapped  phase  and  results  from 
the  direct  and  indirect  approach  are  illustrated  in  Fig.  4.3.  A  linear 
phase  component  is  incorporated  so  that  all  three  impulse  response  esti¬ 
mates  are  causal  and  nonzero  at  n*0.  Note  that  a  jump  of  approximately 
4 it  occurs  near  the  zero  of  H(z)  at  2000  Hz.  This  jump  indicates  that  the 
noise  component  caused  two  encirclements  of  the  origin  by  E(w),  as  depicted 
in  Fig.  4.1b. 

4.2  Phase  Unwrapping  by  Phase-Only  Signal  Reconstruction 

When  H(z)  is  rational,  we  saw  in  section  2.3.2  that  the  number  of  pos¬ 
sible  magnitude  functions  is  finite  when  the  phase  is  fixed.  To  specify 
a  unique  magnitude  additional  constraints  must  be  imposed.  One  set  of  con¬ 
straints  of  particular  interest  is  causality  and  phase  continuity.  These 
constraints  are  compatible  with  only  a  minimum  phase  sequence  and  are  easily 
incorporated  within  a  linear  iterative  procedure  for  magnitude  retrieval. 

A  consequence  of  this  iterative  algorithm  and  the  ultimate  objective  of 


Fig.  4.3  (a)  Unwrapoed  phase  (with  a  linear  phase  component) 
of  H(o))  from  example  4.2 


(b)  Estimate  of  (a)  in  the  presence  of  noise  by  the 
indirect  approach 

(c)  Estimate  of  (a)  in  the  presence  of  noise  by  the 
direct  approach 


I 
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this  section  is  a  new  means  for  computing  the  unwrapped  phase  without 
modulo  2ir  considerations.  Another  consequence  of  this  iterative  proced¬ 
ure,  equally  important,  is  a  new  technique  for  obtaining  the  magnitude 
from  the  phase  of  the  Fourier  transform  of  a  minimum  phase  sequence.  The  ^ 
Hilbert  transform,  the  conventional  approach,  requires  an  unwrapped  phase. 
Our  procedure,  on  the  other  hand,  requires  only  the  principal  value  of  the 
phase,  and  thus  avoids  the  necessity  of  phase  unwrapping. 

Other  constraints  may  also  be  incorporated  within  the  iterative  algor¬ 
ithm.  In  parallel  with  this  thesis,  Hayes  et  al  have  demonstrated  that 
under  a  set  of  rather  loose  conditions  on  H(z)  the  magnitude  function  of 
a  finite  length  sequence  is  unambiguously  determined  (within  a  scale 
factor)  by  the  phase.  As  we  shall  see,  this  finite  length  constraint  can 
therefore  be  imposed  within  the  Iterative  procedure  to  retrieve  the  magni¬ 
tude.  In  addition,  it  is  useful  in  developing  a  OFT  realization  of  the 
iteration. 

4.2.1  Magnitude  Retrieval  From  Phase  With  Constraints 

In  this  section,  we  consider  constraints  for  recovery  of  the  magnitude 
from  the  phase  of  a  Fourier  transform. 

4.2.1. 1  Causality  and  Phase  Continuity 

Consider  a  sequence  h(n)  which  is  causal  and  of  arbitrary  duration. 
Furthermore,  suppose  the  unwrapped  phase  of  its  Fourier  transform  is  con¬ 
tinuous  and  odd,  and  thus  Is  without  a  linear  phase  contribution.  We  assume 
that  h(n)  contains  a  maximum  phase  component  and  then  show  that  this 
assumption  leads  to  a  contradiction. 

h(n)  can  be  expressed  by 


h(">  ’  *  Wn,*hmax(n) 


(4.21) 
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where  hmin(n)  and  hmax(n)  are  normalized  so  that  hm1n(0)»h  (0)«1 . 

Consequently,  the  z-transform  of  h(n)  is  given  by 

H<z)  ’  A  Wz)  Wz) 

'  A(1  + 

x  11  +  nil  lW">z"’>  <4-22> 

Since  the  region  of  convergence  of  Hmax(z)  must  be  the  interior  of  a 
circle  of  finite  radius[22],  the  region  of  convergence  of  H(z)  cannot 
include  z=«.  However,  since  h(n)  is  assumed  causal  and  stable,  from  (4.22) 

lim  H(z)  *  A  (4.23) 

Z-WO 

and  so  a  contradiction  arises,  and  h(n)  must  be  minimum  phase. 

Therefore,  under  our  constraints  of  causality  and  phase  continuity, 
the  magnitude  of  H(u)  is  unambiguously  specified  (within  a  scale  factor)  by 
the  phase. 

4. 2. 1.2  The  Finite  Length  Constraint 

Let  us  suppose  that  H(z)  contains  no  conjugate  reciprocal  zeros.  Under 
this  condition  the  magnitude  of  a  Fourier  transform  is  uniquely  determined 
(within  a  scale  factor)  by  the  phase  when  h(n)  is  constrained  to  be  of 
finite  durat1on[10].  Equivalently,  suppose  h(n)  is  causal  and  of  finite 
duration,  and  H(z)  contains  no  conjugate  reciprocal  zeros.  The  phase  of 
H(w)  is  then  sufficient  to  recover  a  mixed  phase  h(n),  and  the  correspond¬ 
ing  unwrapped  phase  need  not  be  continuous,  i.e.  eh(u>)  may  include  a  linear 
phase  component. 

This  property  can  be  argued  heuristically  from  (2.28).  The  z-transform 
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of  a  finite  length  sequence  contains  only  zeros.  Therefore,  creating  a 
magnitude  function  different  from  the  original  and  maintaining  the  phase 
requires  from  (2.28)  reflecting  a  zero  to  a  conjugate  reciprocal  pole. 

But  the  presence  of  a  pole  implies  an  infinitely  long  sequence,  and  thus 
only  a  zero  geometry  is  allowed.  This  geometry  can  be  shown  to  be  unique 
as  in  [10]. 

Furthermore,  under  these  same  conditions  when  the  sequence  is  M 
points  in  duration,  M  samples  of  the  phase  function  in  the  frequency  inter 
val  [0,ir]  are  sufficient  to  uniquely  determine  the  sequence  within  a  scale 
factor[10]. 

4.2.2  An  Iterative  Procedure  to  Recover  Magnitude  From  Phase 

Solutions  to  both  the  above  magnitude  retrieval  problems  can  be  given 
in  closed  form.  In  the  former  case  the  Hilbert  transform  generates  the 
magnitude  from  the  unwrapped  phase.  In  the  later  case,  the  magnitude  can 
be  obtained  indirectly  through  solution  of  a  set  of  linear  equations[10]. 

In  this  section,  we  discuss  an  alternative  method  of  solution  which 
invokes  a  linear  iterative  algorithm. 

4. 2. 2.1  Theory 

The  iterative  algorithm  will  now  be  described  and  a  proof  is  given 
showing  that  a  defined  error  must  strictly  decrease  as  the  algorithm  iter¬ 
ates.  We  also  prove  that  the  sequence  of  functions  h(c(n)  derived  from  the 
iteration  converge,  to  h(n)  for  nel  where  I  is  the  set  of  integers  for 
which  h(n)  is  known  a  priori.  The  algorithm  is  illustrated  in  Fig.  4.4. 

We  begin  with  an  initial  guess  Mq(u>)  of  the  desired  magnitude  and  take  the 
inverse  Fourier  transform  of  Mg(u)expj[e(u))l  where  e(u)  is  the  known  phase. 
This  step  yields  hg(n),  the  initial  estimate  of  h(n).  Next,  the  known 


Fig.  4.4  Iterative  algorithm  to  recover 
magnitude  from  phase 
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values  of  h(n)  for  nel  are  incorporated  in  the  initial  estimate  (n ) 
to  obtain  h -j (n ) .  The  phase  of  the  Fourier  transform  of  h-j(n)  is  then  re¬ 
placed  by  the  given  phase  and  the  procedure  Is  repeated.  The  steps  Involved 
in  one  Iteration  are  sunmarlzed  below,  hk(n),  0k(u>),  and  Mk(<d)  are  the 
signal,  phase,  and  magnitude  estimates,  respectively  on  the  kth  iteration, 
and  hk+1(n)  is  defined  in  (3.10). 


Iteration  to  Recover  Magnitude  From  Phase 

(i)  Inverse  transform  Mk(u>)  exp[je(u>)]:  (n) 

(ii)  Replace  hk(n)  with  h(n)  for  nel:  hk+i (n ) 

(iii)  Forward  transform  hk+1(n):  Mk+1 (u)exp[j0k+1 (u)] 

(iv)  Replace  0k+1(uj)  with  0(a)):  Mk+1  U)exp[j0(u>)] 

(v)  Repeat 


We  shall  first  show  in  two  steps  that  the  mean  squared  error  between 
H ( u) )  and  HkU)  strictly  decreases  on  each  iteration. 

Error  Reduction  in  the  Time  Domain 

The  mean  squared  error  on  the  kth  iteration  from  Parseval's  Theorem 
can  be  written  as 

IT 

Ek  *  75—  |  | H (o>)  -  Hk(w)|^  dm 

— IT 

"  J  lh(n)  "  hk(n)!2  (4.24) 

3  l  |h(n)  -  h.  (n) [ 2  +  £  |h(n)  -  h.  (n ) | 2 

n^I  nel 
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From  (3.10),  we  have 


|h(n)  -  hk(n)|  *  |h(n)  -  hk+1(n)|  .  "M  (4.25a) 


|h(n)  -  hk(n) j  >  jh(n)  -  hk+1(n)|  =  0  ,  nel  (4.25b) 

Sunning  (4.25a)  and  (4.25b)  over  all  n,  we  obtain 

Ek  -  I  |h(n)  -  hk(n)|2 
K  n  K 


>.  I  |h(n)  -  hk+, (n) 
n  K  1 


(4.26) 


Error  Reduction  in  the  Frequency  Domain 

From  Parseval's  Theorem,  we  write  (4.26)  as 


fit 

Ek  —  5tT  ■  -  ^k+-j(“)l^dw 


(4.27) 


/ir 

|M(u 


2 

)exp[je(u>)]  -  Mk+1((D)exp[j8k+1(o))]|  du 


With  the  triangle  inequality  for  vector  differences,  we  have  at  each  fre¬ 


quency  uj: 


|M(u)exp[je(uj)]  -  Hk+1(«)exp[jek+1(«)]| 

>  |M(u))exp[je(u)]|  -  |Mk+1(o)expCjek+1(a))]| 

=  |M(oj)  -  Mk+^(a))|  (4.28) 

Therefore,  from  (4.27),  (4.28),  and  the  identity  !exp[je(<u)]|2  *  1  : 


^  r 


|M(u)  -  Mk+i  (a>)  |  du> 
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3  ^  \  iM(u>)exp[je(<ii)]  -  (<j))exp[j0(tu)]|^  doi 

-7T 

3  .1. 1 H(<o)  -  H|(+^(uj)|^  du 

_ _ _  -IT 

3  Ek+1  (4.29) 

Therefore,  Ek  Is  at  least  monotone  decreasing  upon  each  Iteration.  Since 
Ek>0,  it  has  a  lower  bound  of  zero.  Therefore, Ek  must  converge  to  a  unique 
limit[30]. 

Let  us  now  suppose  that  on  the  kth  Iteration  hk(n)fh(n)  for  nel. 

From  (4.25b)  and  (4.26)  the  mean  squared  error  must  decrease  on  the  next 
iteration.  Therefore,  the  error  strictly  decreases  unless  convergence  to 
the  limit  of  Ek  has  been  reached. 

Proof  of  Convergence  over  I 

Since  Ek  is  convergent  and  real,  Ek  is  a  Cauchy  sequence[30].  There¬ 
fore,  given  any  real  number  e>0,  there  is  a  positive  integer  N  such  that 

2 

|Ek-EJ<e  whenever  k,*>N,  In  particular,  let  t*k+l  so  that  for  any  e  >0, 
we  can  find  an  N  so  that  whenever  k>N 

IE,,  -  Ek+,|  -  Ek  -  Ek+,  <  c2  (4.30) 

Then,  from  (4.26)  through  (4.29),  and  (4.30) 

l  |h(n)  -  hk(n)|2  -  l  |h(n)  -  hk+1(n)|2 
n  n 

3  l  i h(n)  -  hk(n)|2 
nel 

-  _  2 
<  Ek  "  Ek+1  <  e 


(4.31) 
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Therefore, 

| h (n )  -  hk(n)|2  <  e2  nel  (4.32a) 

and 

|h(n)  -  hk(n)|  <  e  nel  (4.32b) 

Finally,  (4.32b)  Implies  that 

11m  h.  (n)  *  h(n)  nel  (4.33) 

k-*»  K 

From  (4.33),  and  since  each  sequence  hk(n )  has  the  known  phase,  9h(uj), 

a  reasonable  conjecture,  under  the  causality  and  finite  length  constraints 

of  section  4.2.1  where  h(n)»0  for  nel,  is  that 

11m  M.  (ui)  =  aM(oi)  (4.34) 

k-*»  K 

where  a  is  a  constant.  We  shall  see  in  section  4. 2. 2. 4  that  in  practice 
the  limit  given  by  (4.34)  is  in  fact  approached. 

4. 2.2.2  The  DFT  Realization 

We  saw  in  section  4. 2. 1.2  that  M  points  of  the  phase  in  the  frequency 
interval  [0,*]  are  sufficient  to  characterize  an  M  point  sequence  within  a 
scale  factor.  Therefore,  we  are  justified  in  formulating  the  iterative 
algorithm  and  proofs  of  error  reduction  and  convergence  of  section  4.2.1 
for  a  DFT  realization  with  DFT  length  2M-1  or  greater  in  a  manner  similar 
to  that  of  section  3.2. 

4. 2. 2.3  Examples 

The  following  two  examples  illustrate  the  iterative  algorithm  with  the 
constraints  of  section  4.2.1.  In  both  cases  the  initial  magnitude  estimate 
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Mq(u»)  Is  set  to  unity. 

Example  4.3 

Let  h(n)  be  the  sequence  as  defined  in  example  3.1  whose  z- transform 
contains  two  minimum  phase  complex  pole  pairs  and  one  maximum  phase  zero 
pair.  Its  minimum  phase  counterpart  hmp(n)  with  the  continuous  phase  of 
H(oi)  Is  of  Infinite  duration  and  so  within  a  DFT  realization  must  be 
truncated.  We  assume  that  h(n)»0  for  n>256  and  use  a  512-point  DFT.  Trunca¬ 
tion  of  hmp(n)  does  not  noticeably  alter  the  phase  of  H(u>).  Nevertheless, 
the  phase  of  the  Fourier  transform  of  the  truncated  hmp(n)  Is  incorporated 
within  the  iteration. 

The  sequence  of  functions  Mk(o>)  and  e^U)  are  depicted  in  Fig.  4.5 
for  5,  15,  25  and  100  iterations  superimposed  on  the  original  phase,  and 
magnitude  functions  of  hmp(n).  Mk(u>)  and  ek(u>)  are  indistinguishable  from 
the  magnitude  and  phase  of  the  original  minimun  phase  function  after  100 
Iterations.  Note  that  the  maximum  phase  zeros  of  the  original  z- transform 
now  appear  as  conjugate  reciprocal  poles  as  we  would  predict  from  (2.28). 
The  negative  offset  in  the  logmagnitude  estimate  occurs  since  a  in  (4.34) 
is  less  than  unity. 

Example  4.4 

Let  h(n)  be  the  mixed  phase  causal  8-point  sequence  of  example  3.2. 

A  512-point  DFT  Is  used  and  the  linear  phase  component  retained.  The  se¬ 
quence  of  functions  Mk(u>)  and  ek(u)  are  depicted  In  Fig.  4.6  superimposed 
on  the  original  magnitude  and  phase  of  H(u>)  for  1,  5,  15,  and  100  itera¬ 
tions. 


ORIGINAL 


i 

i 
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Fig.  4.5  (a)  Convergence  of  log[Mk(u>)] 

in  example  4.3 
(b)  Convergence  of  0k(o»)  in 
example  4.3  * 
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Fig.  4,6  (a)  Convergence  of  log[M.  (w)] 

In  example  4.4 

(b)  Convergence  of  e.  (u»)  in 
example  4.4 
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4.2.3  A  New  Phase  Unwrapping  Algorithm 

We  are  now  in  a  position  to  state  a  new  phase  "unwrapping"  algorithm 
without  modulo  2ir  considerations.  Suppose  we  are  given  the  principal  value 
of  the  phase  of  H(u>)  and  that  the  corresponding  unwrapped  phase  is  contin¬ 
uous.  Our  new  unwrapping  algorithm  is  outlined  below.  Mmp(ui)  in  step  (i) 
denotes  the  magnitude  of  the  Fourier  transform  of  a  minimum  phase  sequence 
derived  from  our  iterative  procedure  under  a  causality  constraint.  Step 

(i)  yields  the  same  magnitude  (within  a  scale  factor)  that  would  be  obtained 
by  applying  the  Hilbert  transform  to  the  unwrapped  phase,  but  bypasses  the 
need  of  phase  unwrapping. 


Phase  Unwrapping  Algorithm 

(i)  Apply  the  iterative  algorithm  of  section  4. 2. 2.1  with  a  causality 
constraint:  Mmp(<*>) 

(ii)  Compute  the  logarithm  of  the  resulting  magnitude  function  from 
step  (i):  log[Mmp(aj)] 

(iii)  Hilbert  transform  the  logmagnitude  function  derived  in  step  (ii) 
to  obtain  the  desired  unwrapped  phase:  e^U) 


There  are  two  major  considerations  in  the  use  of  this  algorithm. 

First  from  our  discussion  in  section  2.3.2,  we  see  that,  in  general,  the 
minimum  phase  solution  derived  from  the  Iteration  is  of  infinite  extent 
regardless  of  whether  the  original  sequence  h(n)  is  of  finite  duration. 
Therefore,  a  possible  problem  with  aliasing  arises.  The  DFT  length  must  be 
sufficiently  large  so  that  the  minimum  phase  counterpart  of  h(n),  hmp(n) 
has  essentially  decayed  to  zero.  In  particular,  when  h  (n)=0  for  n>N, 
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the  OFT  length  from  the  results  of  section  4,2.1  should  be  at  least  2N-1. 

The  second  consideration  is  the  linear  phase  component  of  H(«i).  The 
presence  of  this  term  represents  a  potential  drawback  to  the  algorithm  since 
a  priori  knowledge  of  such  a  component  is  often  difficult  to  obtain. 
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CHAPTER  5 

ESTIMATION  OF  THE  UNWRAPPED  PHASE  FROM  HARMONIC  SAMPLES 

In  section  4.1.2  we  described  a  direct  approach  to  phase  estimation 
through  homomorphic  deconvolution.  This  approach  relies  on  the  phase  of 
the  quasi -periodic  waveform.  The  magnitude  estimate  is  found  separately 
and  is  not  used  in  deriving  the  phase  estimate.  The  procedure  yields  an 
estimate  which  is  band! imi ted  and  is  an  accurate  representation  of  the 
desired  unwrapped  phase. 

An  implicit  assumption  in  this  direct  approach  is  that  the  pulse 
trainp(n)  of  (2.6)  is  not  exactly  periodic.  When  p(n)  is  periodic  with 
period  P,  P(u>)  is  impulsive  and  periodic  with  period  2tt/P.  Under  this 
condition,  P(u>)  samples  the  desired  system  function  H(oi)  at  harmonically 
related  frequencies,  3  2irk/P,  and  corresponds  to  a  periodic  x(n).  Since 
X(a>)  is  zero  in  frequency  bands,  the  complex  cepstrum  does  not  exist  [33], 
and  thus  our  direct  approach  cannot  be  applied. 

Nevertheless,  this  situation  is  useful  in  representing  voiced  speech 
over  a  short  duration  in  time.  Typically,  within  a  voiced  speech  segment, 
the  vocal  tract  and  vocal  cord  characteristics  are  slowly  varying.  There¬ 
fore,  voiced  speech  over  a  short  duration  (e.g.  20  msec.)  can  be  modeled 
as  a  segment  of  an  infinitely  long  periodic  waveform,  x(n)  3  w(n)  x(n), 
where  w(n)  is  a  finite  length,  unity  amplitude  winaow  over  which  the  vocal 
tract  and  vocal  cords  are  time- invariant. 

Certainly  with  only  samples  of  H(u>) ,  without  additional  constraints, 
there  exist  an  infinite  number  of  choices  for  H(u)  and  likewise  for  the 
phase  of  H(cj).  The  number  of  samples  generally  may  not  be  sufficient  to 


91 


uniquely  specify  the  parameters  of  a  rational  model  through,  for  example, 
techniques  similar  to  those  in  sections  3.2  and  4.2.  Alternatively,  in 
this  chapter  we  take  a  direct  nonparametri c  approach  to  resolving  this 
ambi  gui  ty . 

Specifically,  we  shall  view  the  problem  of  phase  estimation  as  a 
problem  of  polynomial  interpolation  over  harmonic  samples  of  the  real  and 
imaginary  components  of  H(oj),  or  when  only  phase  samples  are  given,  of 
the  all-pass  function  H((jj)/|H(w)  |.  In  particular  we  apply  linear  inter¬ 
polation  over  two  successive  samples.  Conditions  on  9u(w)  are  derived 
under  which  its  corresponding  unwrapped  phase  at  harmonic  samples  is  pre¬ 
served  by  this  simple  procedure  -  a  situation  denoted  as  "phase  tracking". 
Phase  tracking  preserves  the  slowly  varying  component,  i.e.  the  envelope, 
of  the  unwrapped  phase  of  H(ui).  With  an  appropriate  bandlimited  con¬ 
straint  the  entire  phase  function  may  then  be  recovered. 

Linear  interpolation  is  also  useful  in  understanding  the  properties 
of  an  arbitrary  window  w(n)  and  its  relation  to  the  unwrapped  phase  of 
the  windowed  waveform  x(n)  =  w(n)  x(n).  Windowing  can  be  viewed  as  an 
interpolation  procedure  in  the  frequency  domain,  and,  in  fact,  is  itself 
a  means  of  phase  estimation.  With  this  viewpoint,  constraints  on  the 
window  duration  and  alignment  (i.e.  positioning  with  respect  to  x(n)) 
are  imposed  for  phase  tracking. 

In  the  final  section  of  this  chapter,  with  appropriate  a  priori 
information,  a  method  is  derived  for  phase  tracking  by  windowing  without 
the  need  of  alignment. 

5.1  Techniques  of  Phase  Tracking 

Consider  a  periodic  train  of  pulses  p(n)  with  spacing  P.  With 
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p(0)  f  0,  p(n)  in  the  frequency  domain  is  given  by 

P(<*>)  =  p-  I  <S(u-CL)k)  (5.1 ) 

where  =  2ik/P  and  S(ui)  is  the  dirac  delta  function. 

From  section  2.2.3,  we  can  write  the  unwrapped  phase*  of  H(id)  as 

9^(cj)  =  S^(oj)  +  nou  (je[0,7r]  (5.2) 

where  e^w)  is  continuous  and  nQoj  is  a  linear  phase  component.  For  con¬ 
venience,  we  express  the  Fourier  transform  of  h(n)  as  a  function  of  two 
variables: 


H(cj,no)  =  H(u)  exp [jnQo)]  (5.3) 

where 

H(u)  *  | H(ui,  nQ)  |  exp[j0|-(ta))]  (5.4) 

In  the  frequency  domain  we  also  express  the  Fourier  transform  of  the 
periodic  waveform  x(n)  =  p(n)*h(n)  as  a  function  of  two  variables: 


X(to,nQ)  =  H(u),nQ)P(w) 

Therefore,  from  (5.1)  and  (5.5),  we  have 


(5.5) 


X(w,nQ)  =  ^  l  H(ojk,no)  5(ui-u)k)  (5.6) 

and  so  only  samples  of  H(co,n0)  are  available. 

We  have  generally  assumed  throughout  this  thesis  that  h(n)  is 
causal,  nonzero  at  the  origin,  and  has  a  rational  z-transform.  In 
particular,  when  h(n)  is  a  mixed  phase  sequence,  we  see  from  (3.1)  and 
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Throughout  this  chapter  a  phase  function  e(w)  denotes  an  unwrapped 
phase. 
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(3.4)  that  nQ  equals  the  number  of  zeros  of  H(z)  outside  the  unit  circle. 
However,  given  an  arbitrary  segment  of  a  periodic  wavefbrm  of  the  form 
x(n)  *  h(n)*p(n),  we  must  define  an  origin.  Consequently,  without  addi¬ 
tional  knowledge,  the  position  of  h(n)  relative  to  our  defined  origin  Is 
arbitrary,  and  therefore  nQ  In  (5.2)  Is  also  arbitrary.  Because  x(n)  Is 
periodic,  however,  we  can  restrict  nQ  to  the  range  -P/2  <nfl  <  P/2. 

Since,  conceptually  x(n)  Is  of  Infinite  extent,  we  must  apply  a 
window,  w(n).  Throughout  this  chapter,  we  assume  w(n)  Is  symnetrlc  and 
that  the  time  origin  (l.e.  n»0)  is  set  at  the  center  of  w(n).  Therefore, 
changing  nQ  corresponds  to  sliding  x(n)  under  the  window  w(n). 

In  this  section,  we  consider  two  methods  for  phase  tracking,  that  is, 
for  retrieving  the  unwrapped  phase  of  H(w,n0)  at  harmonic  samples.  The 
first  method  of  linear  interpolation  can  be  viewed  as  a  special  case  of 
the  second  method  of  time-domain  windowing.  These  techniques  are  also 
applicable  to  samples  of  the  all-pass  function  H(u,n0)/|H(u,n0) | ,  and 
so  do  not  rely  necessarily  on  magnitude  information. 

5.1.1  Linear  Interpolation  in  the  Frequency  Domain 
One  approach  to  estimating  H(u,nQ),  and  thus  the  phase  %{u)  Is  to 
fit  a  polynomial  of  order  M  to  M  given  samples  of  H(u»k ,nQ )  over  a 
specified  Interval.  For  example,  one  possibility  Is  to  fit  a  first  order 
polynomial  over  two  successive  samples,  i.e.  perform  linear  Interpolation 
on  H^.n^  and  H(ojk+1,n£)).  We  choose  this  particular  scheme  for  three 
reasons  which  will  become  clear  in  the  sequel: 

(1)  Linear  interpolation  illustrates  the  fundamental  problems 
in  preserving  the  unwrapped  phase  of  H(u,n0)  at  harmonics  by 
time-domain  windowing  x(n). 


(II)  The  unwrapped  phase  at  harmonics  derived  by  linear  Inter¬ 
polation  equals  the  unwrapped  phase  found  by  a  running  sum  on 
the  principal  value  of  the  difference  between  two  successive 
samples  of  eh(<u). 

(III)  Phase  estimation  by  linear  interpolation  Is  compatible 
with  the  spectral  envelope  speech  analysis-synthesis  system  given 
In  Chapter  6. 

A  linear  Interpolation  procedure  Is  equivalent  to  convolving  the  real 
and  Imaginary  components  of  X(w,nQ)  with  a  triangular  function  given  by 


W(u»)  • 


"  1  -  |oj/(2tt/P)  | 

< 

0 


H  <  2-n/P 
M  >  2ir/P 


(5.7) 


The  estimate  of  H(u,n0),  H(w,n0)  Is  therefore  expressed  as 

H(u,n0)  -  X(u,nQ)*W(w)  -  7  H(<^.n0)  W(cj-u»|()  (5.8) 


HJ.(w|c*no)  Wtw-u^)  +  j  jp- 1  Hj  (oi|{»no)W(a»-u»|() 


In  the  time  domain,  (5.8)  Is  equivalent  to  multiplying  (l.e.  windowing) 
x(n)  with  a  sequence  given  by 


w(n)  ■  [s1nfern/P)/im]*  (5.9) 

A 

We  say  that  the  unwrapped  phase  of  H(u,n0)  Is  tracked  by  H(u,n0) 
at  w» my  when  the  unwrapped  phase  of  H(u>|(*n0)  equals  the  unwrapped  phase 
of  H(o)jc, nQ): 

8^)  »  (5.10) 


Alternatively,  we  refer  to  (5.10)  as  the  preservation  of  the  unwrapped 
phase  envelope. 
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From  (5.2),  the  phase  difference  A0h(wj()  between  two  successive  samples 
of  0h(<i>)  Is  written  as 

“hK*  *  w  -  9h(“k-i> 

■  "SH*  •  9h(Vl>  +  "o  (“k  •  Vl> 

-  Aflfj  (U|c)  ♦  ZtnJ?  (5.11) 

We  shall  show  that  a  necessary  and  sufficient  condition  for  H(<i>,nQ)  to  track 
the  unwrapped  phase  of  H(u,n0)  at  u-uk  Is  given  by 


0h^“k^  "  flh^“k^  \  W  end  on1y  ^ 

(^(0^)1  »  ♦  2*nQ/P\  <  *\ 

(5.12) 

To  prove  this  condition  we  first  note  from  (5.8)  that  a  linear  trajec¬ 
tory  Is  followed  by  the  real  and  Imaginary  components  of  H(wtnQ)  over  each 
interval  C<*»|c»“|c^3*  Th*  polar  plot  of  H(«,nQ)  for  therefore 

also  follows  a  linear  trajectory  as  depicted  by  trajectories  T1  and  T2  In 
Figs.  5.1a  and  5.1b.  Furthermore,  we  see  from  Fig.  5.1a  that  when 
|A6h(tt|()|<ir,  the  polar  plot  of  H(«,n0)  Is  of  type  Tj,  and  therefore 
A9h^“k^“&flhSc^*  00  ***  other  hand,  when  |a8||(u|C)|>**  the  polar  plot  of 

A 

H(u»,n0)  Is  of  type  T2,  and  therefore  A6|;(«k)«Aeh(»k)+2irMk  where  Is  an 
Integer. 

e£(uk)  can  now  be  expressed  by  the  following  running  sum: 

9s<«k>  *  i0  «(;<“,) 


(b) 

Fig.  5.1  («)  Polar  plot  of  H(«,nQ)  foraeh  <“fc)  •  * 
(b)  Polar  plot  of  H(«,n0)  forae^U^)  >  * 


■  W  +  92,(“k> 


(S.13) 


where  e^ui^)  termed  the  "2*  error  accumulator  function"  Is  given  by 

92,<“k>  *  2”  X  Mt  (5-U) 

Thus  If  | Aeh(u>k) J  <  ir  1^,  from  above  M^O  and  so  from  (5.14) 

Therefore,  from  (5.13)  ®h^“k^“6h^“k^  V  0n  the  other  hand»  1f 
(Ae^toj^)  I >ir  for  some  particular  k,  M^O,  and  from  (5.14)  ®2ir^“k^° 

Finally,  from  (5.13)  *jj(«|()^*||(«|t)  V  and  our  Pr°P°s1tion  is  P»^ven. 

Since,  In  general,  40h(“|{W®h(“|C)+2irM|c  and  from  Fig.  5.1,  l&e^u^)!^, 
we  conclude  from  (2.12b)  that  Ae^u^-PVtae^)]  where  PV  denotes  "prin¬ 
cipal  value  of".  From  5.13  we  then  have 

k 

e^(tuk)  -  PVCA8h(«£)]  (5.15) 

We  see  from  condition  (5.12)  that  the  linear  phase  component  nQu>  plays 
a  major  role  in  phase  tracking.  That  Is,  the  position  of  x(n)  relative  to 
the  window's  center  Is  significant  In  preserving  the  unwrapped  phase  at  har¬ 
monics.  In  Chapter  6,  we  shall  find  that  the  condition  |A0^(uk)|<v  Is 
"almost  always"  true  for  speechllke  harmonics.  Under  this  constraint.  It 
directly  follows  from  (5.12)  that  a  necessary  and  sufficient  condition  for 
phase  tracking  can  be  stated  as 
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^(<\)  *  6h^mv)  **  only  If  nQ  falls  within  the 
range  given  by 

(5.16) 

-  P/2  -  Pm1nAejj(u»k)/2ir  <  nQ  <  P/2  -  PmaxAe^u^)^ 


where  min  and  max  denote  "minimum  or  maximum  value  of"»  respectively. 

When  nQ  falls  within  the  bounds  of  (5.16)  we  say  that  x(n)  Is 
"aligned"  with  our  defined  origin.  If  nQ/0,  but  satisfies  (5.16),  the 
unwrapped  phase  of  0^(0^)  Is  tracked  and  Includes  the  linear  phase  component 
n0uk<  We  can  however  remove  this  term  to  obtain  e^U^)  since  from  (5.2) 


* 


9h(uk1/* 


v* 


(5.17) 


We  thus  far  have  investigated  the  behavior  of  the  unwrapped  phase  of 
H(bi,n0)  at  harmonics.  We  now  shall  demonstrate  that  e^U)  follows  a  mono¬ 
tonic  trajectory  across  each  Interval  [u^  ,wk+ j ] .  To  see  this  we  consider 
the  phase  derivative,  6jj(w)  weighted  by  the  squared  magnitude  of  H(u,nQ). 
From  (2.14),  )H(«,n0) |*8j|(w)  Is  expressed  by 

|H(«,n0)pQ^(<ij)  ■  HjCM^w-w^)  +  Hr(u»|c ,nQ) 3 

(5.18) 

-  MrCM1(u-uk)  +  Hj(«k,n0)],  (acCujj^^] 

1 - 

We  assune  P  Is  even  so  that  H(w,n0)  Is  sampled  at  «■*. 
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* 

where  Mr  and  Mj  are  slopes  of  the  real  and  Imaginary  components  of  H(u»nQ) 

In  the  Interval  Expanding  (5.18),  we  obtain 

|H(w,no}|20j;(a>}  -  Hf(W|(,n0)M1  -  H1(«|(,n0)Mr 

*  c  ,  <i»e[<«i|c,Mjc+1]  (5.19) 

where  c  Is  a  constant.  Since  |H(«,n0)|2  >  0  ,  $£(<■»)  from  (5.19),  Is 
either  always  positive,  negative,  or  zero  In  the  Interval  [<iik,wk+1]. 

Because  e^(u)  Is  given  by  Integration  of  the  phase  derivative,  e^(u)  Is 
therefore  either  monotone  Increasing  or  decreasing  In  our  specified  Inter¬ 
val,  and  we  have  proven  our  proposition. 

Throughout  this  section  we  have  assumed  knowledge  of  samples  of  H(o,n0) 
Me  may  however  apply  linear  Interpolation  directly  to  samples  of  the 
all-pass  function  H(«)/fH(«i)|  and  thus  discard  the  magnitude  Information 
altogether  In  our  procedure.  Clearly,  the  unwrapped  phase  at  samples  de¬ 
rived  from  both  schemes  Is  Identical.  The  difference  In  the  estimate  lies 
In  the  precise  path  of  the  monotonic  trajectory  from  sample  to  sample. 

5. 1.1.1  Examples 

Figure  5.2a  depicts  an  unwrapped  phase  function  e^U)  »  e^U)  +  nQu 
of  ah  all-pass  system  function  H(o>,nQ),  with  no«0,  which  consists  of  factors 
of  the  form  (2,22).  In  the  following  example,  we  consider  estimation  of 
ehU)  by  linear  Interpolation  of  H(u|(,n0)  for  various  sampling  frequencies 
2 ir/P  and  linear  phase  components  n0w.  Any  linear  phase  In  the  estimate 
6£(u)  Is  removed  by  subtracting  U0fl(ir)/ir  to  obtain  an  estimate  of  efjU). 

Example  5.1 

For  a  sampling  frequency  of  2*/P  <  2w/130,  |A6jjUk)|  >  *  at  roughly 


PHASE 


1.5  radians  which  occurs  In  the  decreasing  region  of  depicted  In 
Fig.  5.2a.  Furthermore,  Is  such  that  lae^Cu^)!  <  2ir  so  that 

| A©^(a»|c) |  >  it  for  one  and  only  one  value  of  k  within  this  region.  Figures 
5.2b  and  5.2c  Illustrate  the  estimate  6jj(u>)  (with  any  linear  phase  component 
removed)  for  sampling  frequencies  2ir/130  and  2ir/129,  respectively.  Note 
the  2ir  error  at  approximately  1.5  radians  for  sampling  frequency  2ir/129. 

With  P  *  150,  |aejj{u|()|  <  it  V-^.  Changing  nQ  then  Illustrates  the 
effect  of  the  linear  phase  term  nQu)  In  (5.2).  We  saw  from  (5.16)  that  nQ  is 
constrained  within  a  certain  range  for  phase  tracking.  In  particular,  for 
our  example,  nQ  can  be  shown  to  be  restricted  roughly  to  the  set  n0£[-6,72] 
by  computing  approximate  minimum  and  maximum  Increments  of  ae^U^)  for 
P*150.  Specifically,  consider  the  lower  bound.  Figures  5. 2d  and  5.2e 
Illustrate  the  2ir  error  which  arises  when  nQ  decreases  from  -6  to  -7. 

This  example  Illustrates  the  extreme  sensitivity  of  the  estimate 
e^(ai)  to  changing  either  P  or  nQ  by  as  little  as  one  point, 

5,1.2  Windowing  In  the  Time  Domain 

The  linear  Interpolation  procedure  of  the  previous  section  can  be 
Interpreted  in  the  time  domain  as  multiplication  of  x(n)  by  the  window 
given  In  (5.9).  Consider  now  the  same  window  but  which  Is  not  necessarily 
a  function  of  the  pitch  period: 

wM(n)  -  [slntim/M)/™]2  (5,20) 

where  the  parameter  M  Is  variable.  The  Fourier  transform  of  wM(n), 

WM(a>)  is  given  by  the  triangular  function  of  (5.7)  where  P  is  replaced  by 
M,  Before  proceeding  to  more  general  windows,  we  first  show  that  when 


P  <  M  <  2P,  the  unwrapped  phase  at  harmonics  derived  from  multiplication  by 
wM(n)  Is  Identical  to  that  derived  In  the  frequency  domain  from  the  linear 
Interpolation  scheme  of  the  previous  section. 

A 

Figure  5.3a  Illustrates  the  real  component  of  H(w,n0)  In  (5.8)  where 
successive  weighted  transforms  W^w-o^)  are  partially  overlapping,  l.e. 

P  <  M  <  2P.  In  a  region  of  no  overlap,  because  WM(<o)  Is  real  and  positive 

A 

the  unwrapped  phase  of  H(oi,nQ),  as  depicted  In  Fig.  5.3b,  Is  constant: 

H(w,n0)  *  W^Cw-w^)  | |  <  e/2 

and  so 

*  ^(w^)  Iw-WjJ  <  e/2 

where  e  equals  the  length  of  the  nonoverlapping  region, 
regions  of  no  overlap,  9jj(w)  exhibits  what  we  shall  refer  to  as  "harmonic 
plateaus". 

As  In  section  5.1,1  a  linear  trajectory  Is  followed  by  the  real  and 
Imaginary  components  of  H(u,n0)  In  each  region  of  overlap,  as  depicted  In 
Fig.  5.3a.  Therefore,  a  linear  trajectory  Is  also  followed  by  the  polar 

A  A 

plot  of  H(ai,n0)  in  these  regions.  From  (5.21a),  the  magnitude  of  H(u,n0) 
at  the  endpoints  of  the  overlapping  region  is  scaled  equally.  Consequently, 

A 

the  polar  plot  of  H(w,n0)  is  a  scaled  version  with  the  same  slope  of  that  de¬ 
rived  by  linear  Interpolation  of  harmonics  (i.e.  M*P).  The  unwrapped  phase 
at  each  harmonic  plateau  therefore  must  equal  the  unwrapped  phase  derived 
from  our  procedure  with  M«P.  Thus,  all  properties  given  in  section  5.1.1 
hold  also  for  P  <  M  <  2P. 


(5.21a) 

(5.21b) 
Therefore,  in 


Let  us  now  suppose  that  we  are  given  an  arbitrary  finite  length  sym¬ 
metric  window,  with  Fourier  transform  W(w).  From  (5.8)  we  can  view  window- 
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Ing  Itself  as  a  phase  estimation  procedure.  The  phase  trajectory  between 
harmonics  Is  determined  by  the  malnlobe  and  sldelobe  structure  of  W(u>)  and 
the  degree  to  which  and  WU-w^)  overlap.  If  we  approximate  the 

A 

trajectory  of  H(a»,nQ)  as  linear  in  a  region  where  malnlobes  "significantly" 
overlap,  and  If  we  assume  that  outside  this  region  the  phase  is  constant  as 
in  (5.21b),  then  the  phase  trajectory  follows  that  derived  with  wM(n). 

For  example,  the  Fourier  transform  of  a  Hamming  window  has  a  malnlobe 
which  can  be  approximated  roughly  by  the  triangular  function  WM(w).  Fur¬ 
thermore,  the  duration  of  a  Hamming  window  should  be  less  than  about  four 
pitch  periods.  Otherwise  there  Is  no  region  of  mainlobe  overlap  and  the 
phase  trajectory  Is  determined  primarily  by  sldelobe  structure,  and  weights 
H(“k*no)  and  H^“k+l’no^*  Me  ^ave  f°und  such  approximations  to  be  useful  in 
describing  the  empirical  behavior  of  the  unwrapped  phase  of  a  periodic  wave¬ 
form  multiplied  by  a  Hamming  window. 

5. 1,2.1  Examples 

In  the  following  example  we  consider  the  estimation  of  the  unwrapped 
phase  function  9h(o»)*0^(u)  +  nQu  of  H(u>)  from  example  4.1,  and  which  is  de¬ 
picted  in  Fig.  5,4a  with  nQ=0,  We  apply  a  Hamming  window  and  demonstrate 
the  effect  of  modifying  its  duration  and  and  its  position  with  respect  to 
the  periodic  waveform  x(n)«h(n)  *  p(n)  where  p(n)  has  a  period  of  50  points. 
9~(u>)  is  such  that  | |  <  it, 

Example  5.2 

Figures  5.4b  and  5.4c  depict  the  continuous  unwrapped  phase  estimates 
obtained  with  application  of  a  Hamnlng  window  of  length  2  and  3.9  pitch 
periods,  respectively,  and  where  nQ»0.  Note  that  when  P  equals  two  pitch 
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iSVHd 


Fig.  5.4  (a)  Unwrapped  phase  of  a  system  function  with  two  complex  pole  pairs  and  one 
complex  zero  pair,  (b)  Estimate  of  (a)  In  example  5.2  with  a  100  point  Hanning 
window  where  n  =0,  (c)  Same  as  (b)  with  a  390  point  Hamming  window,  (d)  Error 
function  for  estimate  of  (a)  In  example  5.2  with  a  100  point  Hanning  window 
where  n  =-20. 
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periods,  the  phase  trajectory  Is  generally  monotonlc  between  regions  where 
the  phase  Is  roughly  constant  near  harmonics. 

Figure  5.4d  depicts  the  error  function  e  (w)  given  by 

8e(<*0  ■  -  9hU)  (5.22) 

for  the  unwrapped  phase  derived  with  a  Hamming  window  of  duration  two  pitch 
'  periods,  but  where  nQ  ■  -20.  (S.22)  may  be  viewed  as  a  continuous  counter¬ 

part  to  the  2ir  error  accumulator  function  of  (5.14)  where  0e(“j()  * 

From  Fig.  5.4d  we  deduce  that  u»n0  Is  such  that  lae^U^)  |" | A©*(ai|t) 

+  2irnQ/P|  >  it  In  the  region  of  a  minimum  phase  pole  (denoted  by  P)  and  a 
maximum  phase  zero  (denoted  by  Z).  Consequently,  there  arises  an  error  of 
2ir  In  both  regions. 

5.1.3  Comments  on  Comparing  Magnitude  and  Phase  Estimation  by  Windowing 

We  have  seen  In  sections  5.1.1  and  5.1.2  that  linear  Interpolation  of 
samples  of  a  Fourier  transform  through  windowing  with  w^n)  generates  a 
monotonlc  phase  trajectory  In  the  interval  Such  monotonicity  is 

typical  of  speechlike  phase  between  harmonics.  That  is,  we  do  not  expect 
"ripples"  in  the  phase  function,  but  a  "smooth"  behavior  between  successive 
harmonics. 

Let  us  now  address  magnitude  estimation  by  this  Interpolation  procedure. 
In  particular,  let  us  suppose  that  we  are  given  two  successive  samples  of 
H(iu,n0),  one  to  the  right  and  one  to  the  left  of  a  pole  located  close  to 
the  unit  circle.  Furthermore,  we  assume  nQa0  and  that  |AehUk)| 

■  |A0^(uk)|  =  IT. 

From  Fig.  5.1,  we  see  that  the  magnitude  of  H(w,n0)  generates  a  null 
where  In  fact  there  exists  a  peak  In  H(w,n0)  due  to  the  presence  of  a  pole. 
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Such  "pole  splitting"  Is  therefore  potentially  a  problem  with  this  tech¬ 
nique.  Linear  Interpolation  applied  directly  to  magnitude  samples  may  thus 
be  preferred  In  such  cases. 

5.2  On  The  Problem  of  Alignment 

In  section  5.1.1  we  derived  the  bound  (5,16)  on  n0  to  ensure  phase 
tracking  under  the  constraint  (ae^U^)!  <  v  We  referred  to  this  con¬ 
dition  as  alignment  of  x(n).  A  number  of  examples  were  given  In  the  pre¬ 
vious  sections  Illustrating  the  Importance  of  this  requirement.  Our  main 
purpose  In  this  section  Is  to  demonstrate  that  alignment  Is  an  Inherently 
ambiguous  process  with  the  sole  constraint  that  j A0^(ti»k)  |  <  ir.  Thus  without 
additional  a  priori  knowledge,  an  alignment  guess  must  be  made. 

In  section  5.2.2  a  method  Is  described  for  phase  tracking  without  the 
need  of  alignment,  under  appropriate  constraints  on  the  second  difference 
of  ) . 

5  2,1  Alignment  Antiquity 

With  the  constraint  lae^O^)!  <  ir  V>k,  we  would  like  to  determine  by 
Inspection  of  0^(0^)  whether  nQ  falls  within  the  bounds  of  (5.16)  so  that 
phase  tracking  Is  guaranteed.  From  (5.13)  our  problem  Is  equivalent  to 
detecting  the  presence  of  ^^“k^  wh1ch  *s  zer0  ^  and  only  n0 
(5.16). 

One  situation  which  can  arise  and  lead  to  an  erroneous  conclusion 
about  alignment  Is  the  case  where  e^(ir)«0  which  occurs  when  nQ»0.  However, 
suppose  nQ  »  -  Then  from  (5,13)  the  2ir  error  accumulator  function 

9 2 w ( <*>k ^  canc®1s  the  linear  phase  component  at  w«it.  Therefore,  e^(ir)»0  Is 
not  a  necessary  condition  for  no"0,  nor  for  phase  tracking. 


This  case  Is  Indicative  of  the  problem  of  alignment  ambiguity.  More 
generally,  we  see  from  Fig.  5.1  that  without  additional  a  priori  knowledge, 
we  cannot  determine  whether  the  true  net  phase  trajectory  between  harmonics 


Is  clockwise  or  counterclockwise.  This  ambiguity  arises  since  when  nQu  Is 
sufficiently  large,  |&eh(uk)|  >  it. 

In  conclusion,  we  cannot  detect  the  presence  of  ^(o^)  ^ro®  only 

0h(“k)* 

5.2.2  Phase  Tracking  Without  Alignment 

One  method  of  phase  tracking  which  does  not  depend  on  alignment  capi¬ 
talizes  on  the  constraint  that  the  second  difference  of  e^U^)  does  not 
exceed  it.  nils  constraint  Is  natural  for  speechllke  phase  because  It  re¬ 
moves  the  possibility  of  "large'1  ripples  In  the  unwrapped  phase  function. 
Our  technique  provides  a  way  to  undo  the  effect  of  misalignment  and  relies 
primarily  on  the  fact  that  second  differencing  of  samples  of  the  unwrapped 
phase  effectively  eliminates  the  linear  phase  contribution  at  each  sample. 
With  some  algebraic  manipulation  the  second  difference  of  can 

be  expressed  by 

2 

A  9|j(uij^)  *  A6^(<i>|()  -  )  (5.23) 


A_ec(u|.) 


nowl 5 (“k 


U1 


) 


*  92,(“k>  -  292,(Vl>  *  92^(“tc-2) 

and  so  the  linear  phase  component  Is  removed  except  at  u>*u<|s2ir/P. 
With  the  constraints. 
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2*n_ 

!•*(«,)  ♦-t£I  <  *  (5.24b) 

and,  the  Initial  condition, 

AjjUj)  -  Ae^(a»7 )  «  ejjUj)  (5.25) 

2 

we  can  write  the  principal  value  of  A  e£(«k)  as 


2  2  9 

A  *jj{wk)  ■  PV[AS|j(«k)]  -  A  »f;(<*»k)  ♦  — p-1  «(u»k-%*7 )  (5.26) 

The  running  sum  of  (5.26)  Is  taken  to  obtain 

2im 

AljjUjj)  *  Afl|j(wk)  +  p  (5.27) 

Repeating  for  (5.27),  we  have  the  desired  unwrapped  phase  at  harmonics 

♦'(a^)  *  »  0jj(»k)  *  no"k  (5.28) 

2 

Note  that  PV  [a  e^c*^}]  need  not  be  computed  by  unwrapping  and  differencing, 
but  rather  directly  through 


PV  [A20£(<nk)]  ■  tan"1 


s1n[AZ0j;(wk)] 

COS[AZ0£(wk)] 


(5.29) 


2  2 

since  s1n[A  e^(uk)]  and  cos[a  9£(wk)]  can  be  expressed  In  terms  of 
cos[8|j(uk)]  and  s1n[8£(uk)]. 
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CHAPTER  6 

APPLICATION  OF  PHASE  TRACKING  TO 
SPEECH  ANALYSIS-SYNTHESIS 

In  this  chapter,  we  apply  tha  tachnlquas  of  phasa  tracking  of  tha  pre¬ 
vious  chaptar  In  Introducing  a  alxad  phasa  astlnata  within  two  spaach  anal- 
ysls-synthasls  systeam: 

(I)  Tha  homomorphic  systam  propostd  by  0ppenhe1mC23] 

(II)  Tha  spactral  anvalopa  systam  proposad  by  Paul [26] 

With  a  slmpla  modal  of  tha  phasa  of  tha  vocal  tract  frequency  response,  we 
first  demonstrate  that  speechllke  bandwldths  are  such  that  an  unwrapped 
phasa  Increment  across  successive  harmonics  generally  does  not  exceed  *  , 
l.e. ,  <  v.  Many  of  tha  results  of  tha  previous  chaptar  are  there¬ 

fore  applicable  In  tha  analysis  of  voiced  spaach  whan  modeled  over  a  short 
duration  as  a  segment  of  a  periodic  waveform. 

Tha  concept  of  short-time  homomorphic  analysis  of  a  periodic  waveform 
Is  Introduced  and  Is  shown  to  rely  strongly  on  tha  nature  of  tha  window  In 
tha  time  domain.  A  Hamming  window,  with  a  pitch-adapted  duration,  and 
appropriately  aligned.  Is  applied  to  Improve  the  phase  estimate  derived  from 
homomorphic  deconvolution.  Such  a  windowing  procedure,  results  In  the  addi¬ 
tional  property  that  our  analysis-synthesis  system  Is,  potentially,  an 
Identity  system  with  respect  to  a  periodic  waveform. 

Finally,  linear  Interpolation  of  section  (5.1.1)  Is  directly  applied 
In  deriving  a  phase  estimate  which  Is  Incorporated  within  the  spectral  en¬ 
velope  system.  The  original  scheme  based  on  magnitude  only  relies  on  deter¬ 
mining  harmonic  locations  by  peak  picking,  and  thus  Is  compatible  with  our 
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Interpolation  procedure  for  phase  tracking. 

Informal  listening  tests  Indicate  a  small  but  perceptible  Improvement 
In  "quality"  within  these  systems  when  a  mixed  phase  estimate  replaces  the 
minimum  phase  counterpart. 

6.1  Bandwidth- Pitch  Period  Constraints 

The  maximum  phase  Increment  due  to  a  single  pole  or  zero  approaches  * 
as  either  comes  close  to  the  unit  circle.  The  Increment  will  not  exceed  *, 
so  that  In  these  simple  cases  |ae£(u|()|  <  t.  The  geometry  of  a  more  compli¬ 
cated  pole- zero  pattern,  however,  might  be  such  that  |ae^(«^)|  >  v  .  It  Is 
the  purpose  of  this  section  to  derive  bounds  on  the  bandwidth  of  poles  and 
zeros  of  a  z-transform  for  which  this  constraint  Is  satisfied,  and  to  show 
that  typical  speech  bandwldths  fall  within  these  bounds. 

Fig.  6.1  depicts  two  regions,  A  and  B,  of  an  elemental  phase  function, 
e(u)  which  corresponds  to  a  real  maximum  phase  zero  with  Fourier  transform 
H(u)  »  1  -  aexpCJu],  0<a<1  .  Any  complex  pole  or  zero  contributes  a  phase 
which  Is  a  shifted  and/or  negated  version  of  that  In  Fig.  6.1.  An  arbitrary 
phase  Increment  ae^u^)  can,  therefore,  be  represented  by  a  sum  of  Incre¬ 
ments  derived  from  e(u). 

To  obtain  bounds  on  this  sum,  we  consider  the  phase  derivative  of  e(u): 
©  2  2 

eU)  ■  (a  -  acos«)/(a  -  2acos«  *1)  (6.1) 

The  maximum  of  e(»)  In  region  B  occurs  at  w  *  +  ir  : 

max  e(w)|  ■  ${w)|  ■  a/(a+l)  (6.2) 

'B  '  Ui*+TT 

and  the  minimum  In  region  A  occurs  at  u  ■  0: 
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min  «(«)!  *  9(«) 


a/(a-l) 


(6.3) 


For  simplicity  and  mathematical  tractablllty,  we  express  a  phase  Incre¬ 
ment  over  the  Interval  as  a  sum  of  Increments  from  region  A  and 

from  region  B.  A  minimum  phase  pole  or  maximum  phase  zero  contributes  a 
phase  of  the  form  e{«-«i»0).  A  minimum  phase  zero  contributes  a  phase  of  the 
form  >e(u-u0)  so  that  the  contribution  to  a  phase  Increment  from  region  A 
or  B  may  therefore  be  negated. 

Since  e(b>)  Is  monotonlc  in  region  A  or  B,  from  (6.2)  and  (6.3)  bounds 
on  an  Increment  A9(wk)  In  the  two  regions  are  given  by: 


0  >  > -Stan'1  [r-  > 


(6.4) 


0  <  dd(u|^)  J 


ir/P  <  it 


(6.5) 


Me  consider  two  cases  which  are  clearly  not  exhaustive,  but  Indicative 
of  bandwidth  constraints: 

Case  1  (sparsely  spaced  poles  and  zeros): 


When  the  poles  and  zeros  of  H(z)  are  sparsely  spaced,  we  shall  assume 
that  any  Interval  of  8jj(«)  consists  of  either  the  sum  of  contributions  from 
only  B  regions,  or  the  sum  of  contributions  from  all  B  regions  and  at  most 
one  A  region. 

Mhen  only  6  regions  overlap,  from  (6.5)  !<^(<»|()l  <  *  because  we  assume 
the  nunber  of  poles  and  zeros  to  be  less  than  P.  When  there  exists  one  A 
region,  we  assume  the  sum  of  contributions  from  B  regions  to  be  positive 
and  less  than  *.  The  positivity  constraint  holds  since  for  speechlike 


114 


spectra,  the  number  of  poles  Is  greater  than  the  number  of  zeros.  There* 

fore,  since  0  >  40(0^)  I  >  -  *,  we  have  |40|j(o^)|  <  »  . 

'A  f1 

Case  2  (closely  spaced  poles  and  zeros): 

Consider  the  overlap  of  two  A  regions  (e.g.,  two  closely  spaced  poles  or 

a  closely  spaced  pole  and  maximum  phase  zero)  and  an  arbitrary  number  of  B 

regions.  Since  we  assume  the  net  contribution  of  the  B  regions  Is  positive, 

from  (6.4)  and  (6.5),  we  write  the  following  approximate  constraint 

l“h("k>l<l2  t,n‘1  Irf  losiffljj  <6-6> 

where  IOr/2P  Is  the  "average"  contribution  from  K  B  regions  and 
2  tan"1 

is  the  "average"  contribution  from  two  A  regions,  where  we  have  assumed 
equal  variable  bandwldths. 

With  some  algebraic  manipulation  our  constraint  can  be  expressed  In 
terms  of  half  power  bandwidth  a [5]  by 

a  -  ^Lnja"1 

>  ^-Ln  cos(ir/P)  +  sin(ir/P)  tanC (J(  1+K/2P) J”1  ]  (6.7) 

where  F$  Is  the  A/0  sampling  frequency. 

Table  6.1  gives  values  of  a  for  typical  values  of  P  and  possible  values 
of  K  corresponding  to  four  or  five  real  or  complex  pole  pairs,  and  one  or 
two  real  or  complex  zero  pairs.  Table  6.2  gives  average  measured  bandwldths 
for  the  first  three  formants  of  vowel  utterances  for  three  male  subjects 
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TABLE  6.1 

Computed  half-power  bandwldths  required  for  phase  tracking 


K 

P 

1 

2 

3 

4 

5 

6 

7 

8 

9 

10 

4  msec 

83 

82 

80 

79 

78 

77 

76 

74 

73 

72 

6  msec 

58 

57 

56 

56 

55 

54 

54 

53 

53 

52 

8  msec 

44 

43 

43 

42 

42 

42 

41 

41 

41 

40 

10  msec 

36 

35 

35 

34 

34 

34 

34 

33 

33 

33 

TABLE  6.2 

Mean  half-power  bandwldths  (Bl,  B2,  B3)  In  cycles  per  second 
for  the  first  three  resonances  of  vocal  tract  configurations. 
Three  male  subjects  (SI,  S2,  S3)  used  eight  vowel  conflgura- 
tlons  with  two  glottal  conditions  (after  House  and  Stevens[ll J). 


Band 

Width 

SI 

Open  Glottis 

S2  S3  Mean 

SI 

Closed  Glottis 

S2  S3  Mean 

Bl 

79 

75 

67 

73 

56 

50 

56 

54 

B2 

88 

82 

75 

81 

67 

62 

61 

65 

B3 

97 

98 

83 

94 

76 

66 

64 

70 
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with  open  and  closed  g1ott1s[11].  Since  during  a  vowel  utterance  the  glottis 
Is  partly  open  and  partly  closed,  the  actual  bandwldths  probably  lie  some* 
where  between  those  measured  for  the  open  and  closed  glottis.  It  Is  clear 
that  these  bandwldths  often  fall  within  the  derived  lower  bounds.  It  Is 
conceivable,  however,  that  spectra  of  high-pitched  speakers  with  closely 
spaced  formants  may  yield  a  phase  Increment  outside  of  our  desired  constraint 
(l.e.,  |A6|j(tok)  |<ir). 

6.2  Pre-  and  Post-Al i qnment 

Alignment  of  x(n)  ■  h(n)*p(n)  so  that  nQ  falls  within  the  constraints  of 
(5.16)  will  be  referred  to  as  pre-alignment.  Under  this  constraint 
iAeh(“k)l<ir  so  that  w-indowin9  the  t-ime  domain  or  performing  linear  inter¬ 
polation  in  the  frequency  domain  tracks  the  desired  unwrapped  phase  at  har¬ 
monics,  Since  n0*e^(ir)/ir,  the  linear  phase  term  nQu>  can  be  eliminated  by 
subtraction  —  a  procedure  we  shall  refer  to  as  post-alignment.  This  oper¬ 
ation  generates  an  estimate  of  e^(<i>)  which  equals  the  desired  unwrapped 
phase  at  harmonics  when  (5.16)  holds,  and  is  an  accurate  approximation  to 
the  desired  unwrapped  phase  between  harmonics.  Although  pre-alignment  may 
guarantee  that  (5,16)  holds,  the  value  of  nQ  may  not  be  consistent  from  frame- 
to-frame.  Consequently,  when  nQ  Is  not  removed  "pitch  jitter"  can  arise  at 
frame  boundaries  after  signal  reconstruction,  causing  degradation  in  the 
synthesized  speech.  That  Is,  the  arbitrary  displacement  of  the  vocal  tract 
Impulse  response  estimate  results  In  a  random  change  in  pitch  from  frame- 
to-frame. 

In  practice,  however,  estimation  and  removal  of  nQ  Is  complicated  by 
deviations  from  the  assumed  harmonic  structure  of  X(u»,n0)  In  (5.6)  In  high- 
frequency  regions.  Two  causes  for  such  deviations  are:  (1)  low-pass  filter- 
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ing  the  speech  waveform  before  A/D  conversion,  and  (il)  a  mixed-source 
excitation  at  the  glottis  which  is  harmonic  in  a  low-frequency  region  and 
flat  In  a  high-frequency  reg1on[16].  Therefore,  deriving  nQ  from  (w) 
is  not  recommended. 

One  simple  way  of  representing  high-frequency  degradation  In  the  case 
of  low-pass  filtering  is  through  a  model  of  x(n)  given  by 

x(n)  »  h(n)*g(n)*p(n)  (6.8) 

where  g(n)  is  a  low-pass  filter.  Ideally,  g(n)  introduces  a  linear  phase 
component  within  the  phase  of  H(w).  The  high-frequency  energy  of  g(n)  is, 
however,  quite  low  and  so  in  practice,  the  unwrapped  phase  of  g(n)  beyond 
its  low-frequency  cutoff  Is  erratic  due,  for  example,  to  quantization  noise. 
Consequently,  the  value  of  is  unpredictable  and  may  have  no  relation¬ 

ship  to  the  original  linear  phase  of  H(w).  Subtraction  of  aje^(ir)/Tr,  there¬ 
fore,  arbitrarily  shifts  the  estimate  of  h(n)  and  results  in  "pitch  jitter" 
at  frame  boundaries.  There  exist  then  two  causes  of  "pitch  jitter":  (i) 
inconsistent  pre-alignment,  and  (ii)  erroneous  post-alignment. 

For  the  purpose  of  waveform  reconstruction,  it  is  necessary  to  obtain 
only  relative  delays  between  successive  impulse  response  estimates.  An 
alternative  method  of  post-alignment  invokes  the  cross-correlation  function 
of  two  successive  estimates.  Given  that  two  successive  estimates  hk(n) 

A 

and  hk+-j(n)  are  slowly  varying  except  for  a  delay,  their  cross-correlation 
function  Is  given  roughly  by 

R(n)  =  hk(n)*hk+1(-n) 

=  [h(n)*«(n-nk)*g(n)]*[h(-n)*5(-n-nk+1  )*g(-n)] 

■  Rfj(n)*Rg(n)*s(n-(nk-nk+1 )) 


(6.9) 
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where  (n )  and  R^(n)  are  the  autocorrelation  functions  of  h(n),  our  de¬ 
sired  response,  and  g(n),  respectively,  and  and  nk+^  are  the  delays  in 
h(n)  due  to  the  pre-alignment  process.  Therefore,  the  location  of  the  peak 
in  R(n)  is  an  estimate  of  n^-n^],  and  does  not  rely  on  the  unwrapped  phase 
at  u*ir.  Me  can  then  perform  accurate  post-alignment  by  shifting  the  (k+l)st 
impulse  response  estimate  by  n^-n^  points. 

6.3  Homomorphic  Speech  Analysis-Synthesis 

In  this  section  we  first  review  the  minimum  phase  homomorphic  speech 
analysis-synthesis  system  proposed  by  Oppenheim.  With  the  results  of 
Chapter  5  and  sections  6.1  and  6.2  as  a  foundation,  we  then  develop  a  high- 
quality  homomorphic  system  which  incorporates  a  mixed  phase  estimate. 

6.3.1  Minimum  Phase  Analysis- Synthesis 

The  analyzer  of  the  minimum  phase  system  consists  of  Fourier  trans¬ 
forming  a  short-time  speech  segnent,  computing  the  logmagnitude  of  its 
Fourier  transform,  and  inverse  transforming  to  generate  the  real  cepstrum. 
Pitch  information  is  obtained  with  a  cepstral  pitch  detector[3]  which 
utilizes  cepstral  peak-picking,  and  energy  and  zero-crossing  measurements. 
Three  and  five  point  median  smoothing[3]  was  applied  to  estimate  one  and 
two  point  Isolated  singularities  due  to  pitch  doubling  and  halving,  or 
voiced/ unvoiced  errors.  Some  hand  editing  was  performed  at  voiced/unvoiced 
transitions. 

The  minimum  phase  Impulse  response  estimate  is  derived  by  multiplying 
the  real  cepstrum  with  a  3.2  msec,  low-time  gate  of  the  form  in  (4.15). 

The  result  is  then  transformed,  exponentiated,  and  Inverse  transformed. 

The  waveform  is  synthesized  by  explicitly  convolving  the  impulse  re¬ 
sponse  estimate  and  excitation.  During  voicing,  the  excitation  consists 
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of  a  train  of  unit  inpulses  with  spacing  equal  to  the  pitch  per*<*d.  During 
unvoiced  intervals,  a  noiselike  waveform  of  random  polarity  with  spacing 
1  msec,  is  used.  Linear  interpolation  of  the  pitch  period  and  impulse 
response  was  performed  to  avoid  sudden  changes  in  pitch  and  spectral  infor¬ 
mation.  Such  interpolation  leads  to  enhanced  quality [23, 28].  This  system 
was  simulated  on  a  PDP/11 -55  computer  with  floating  point  arithmetic. 
6db/octave  pre-  and  de-emphasis  was  used,  and  the  input  speech  was  low-pass 
filtered  at  4.8  kHz. 

The  resulting  synthetic  speech  Is  of  high  quality,  natural  sounding,  and 
provides  a  reference  point  for  the  mixed  phase  simulations  to  follow. 

6.3.2  Mixed  Phase  Analysis 

The  first  attempt  to  Introduce  a  mixed  phase  estimate  within  the  homo¬ 
morphic  scheme  through  the  complex  cepstrum,  as  defined  through  (4.1)  and 
(4.2)  resulted  in  large  sensitivity  of  the  phase  estimate  to  the  position 
and  duration  of  the  time-domain  window£31 ,35J,  The  system  generated  synthetic 
speech  characterized  by  a  "hoarse"  quality.  This  sensitivity  Is  clear  in 
light  of  section  5.1. 

The  reason  for  the  hoarseness  is  that  neither  pre-  nor  post-alignment 
was  performed  and  the  window  duration  was  In  general  greater  than  four 
times  the  pitch  period.  Consequently,  both  phase  degradation  and  pitch 
jitter  are  Introduced  within  the  reconstructed  waveform.  Furthermore,  our 
results  illustrate  that  smoothing  the  phase  through  the  complex  cepstrum 
Is  meaningful  only  when  the  envelope  of  the  phase  of  H(a»,nQ)  is  preserved; 
that  Is,  when  phase  tracking  is  guaranteed  by  appropriate  windowing.  The 
requirements  on  the  window,  however,  are  such  that  the  windowed  speech 
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waveform  may  not  follow  a  convolutional  model [31 ,35]  —  a  phenomenon 
which  can  be  easily  seen  from  the  results  of  section  5,1.  Therefore,  pitch 
and  Impulse  response  Information  In  the  complex  cepstrum  are  not  neces¬ 
sarily  additive.  Empirical  results  Indicate  that  for  an  arbitrary  window 
the  breakdown  of  the  additivity  assumption  occurs  severly  with  respect  to 
phase,  but  far  less  so  with  respect  to  magnitude. 

Applying  a  low-time  gate  to  the  complex  cepstrum  may  be  alternatively 
Interpreted  as  a  method  of  bandllmltlng  the  complex  logarithm,  and  thus 
smoothing  the  logmagnltude  and  phase.  When  the  phase  Is  not  tracked  at  har¬ 
monics  as,  for  example,  when  x(n)  does  not  fall  within  the  constraints  of 
(5.16)  or  the  window  Is  excessively  long,  large  erroneous  trends  are  Intro¬ 
duced  Into  the  phase  function  as  depicted  In  Figs.  5.2  and  5.4.  In  the 
former  case,  these  trends  are  due  largely  to  the  2*  error  accumulator  func¬ 
tion,  92ir(«).  It  is  reasonable  to  assume  that  smoothing  e2ir(<i))  can  alter  the 
principal  value  of  the  desired  phase  in  high  energy  regions  or  at  harmonics. 
Therefore,  with  respect  to  either  h(n)  or  a  reconstruction  of  the  original 
waveform  x(n),  signal  distortion  due  to  phase  distortion  may  be  large. 

Such  distortion  has,  in  fact,  been  observed  empirically  in  both  the  estimate 
of  h(n) [28,35]  and  the  reconstructed  waveform. 

6.3.3  Short-Time  Reconstruction 

The  notion  of  bandllmltlng  the  complex  logarithm  is  an  interesting  one 
since  It  Implies  an  Important  alternative  interpretation  of  the  synthesis 
procedure.  As  the  cepstral  window  approaches  unity  the  original  windowed 
waveform  Is  preserved.  Therefore,  when  the  windowed  waveform  has  primarily 
a  low-time  cepstral  composition,  gating  its  complex  cepstrum  maintains  Its 
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basic  characteristics. 

Furthermore,  when  a  Hamming  window  Is  two  pitch  periods  In  extent,  and 
positioned  Identically  with  respect  to  h(n)  In  successive  frames,  the  syn¬ 
thesis  procedure  reconstructs  a  periodic  waveform  exactly  and  thus  acts  as 
an  Identity  system.  This  Is  because  a  Hanning  window  of  two  pitch  periods 
In  length  has  zeros  at  the  spectral  harmonics.  Equivalently,  a  Hamming  win¬ 
dow  when  repeatedly  added  to  Itself  delayed  by  a  pitch  period  results  In 
unity[28].  The  Identity  system,  however,  does  not  rely  on  satisfying  the 
alignment  requirement  of  (5.16)  but  simply  a  consistent  positioning  of  the 
window  from  frame- to-frame.  The  overall  analysis-synthesis  system  Is  de¬ 
picted  In  Fig.  6.2. 

An  example  of  the  output  of  this  system  for  a  real  speech  input  is 
illustrated  In  Fig.  6.3  for  51.2  msec,  of  data  across  three  frame  boundaries 
with  a  20  msec,  frame  rate.  The  cepstral  window  In  this  particular  case  is 
of  duration  .8  P  (where  P  equals  the  pitch  period),  demonstrating  the  scheme's 
potential  as  an  Identity  system.  The  pre-alignment  process  Is  heuristic, 
consisting  of  picking  the  maximum  absolute  value  within  the  second  pitch 
period  from  the  current  frame  number  and  back  tracking  10%  of  the  pitch 
period.  The  synthesis  includes  also  post-alignment  by  cross-correlation 
of  section  6.3. 

6.3.4  Informal  Listening  Tests 

An  analysis-synthesis  system  incorporating  a  mixed  phase  estimate  was 
designed  and  evaluated.  A  Hanning  window  of  two  pitch  periods  In  duration 
was  applied  In  the  analysis.  Our  requirements  on  window  length,  pre-align¬ 
ment  and  post-alignment  for  periodic  segments  reduced  the  characteristic 
hoarseness  of  this  system.  Removal  of  any  one  of  these  requirements  Increased 
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hoarseness.  A  fixed  20  msec,  time  window  was  applied  to  unvoiced  segments, 
and  a  symmetric  rectangular  low-time  gate  with  cutoff  3.2  msec,  was  applied 
to  the  complex  cepstrum. 

When  compared  with  Its  mlnlmun  phase  counterpart,  the  system  with  mixed 
phase  produced  small,  but  audible.  Improvement  In  Uqua11tyn.  The  two  systems 
were  Identical  except  for  the  method  of  Introducing  phase  and  the  elimination 
of  post-alignment  In  the  minimum  phase  version. 

An  Informal  A/B  listening  test  was  performed  where  each  listener  was 
asked  to  choose  the  synthetic  speech  passage  closest  to  the  original.  Ten 
listeners  and  eight  sentences  with  5  male  and  3  female  speakers  were  used. 

The  system  with  mixed  phase  was  judged  roughly  45%  of  the  time  to  be 
closer  In  quality  to  the  original  than  Its  minimum  phase  counterpart,  the 
minimum  phase  version  10%  of  the  time  was  judged  closer,  and  45%  of  the  time 
the  two  were  Indistinguishable.  When  preferred,  the  system  with  phase  was 
often  judged  by  experienced  listeners,  to  reduce  "buzzlness"  of  the  mini¬ 
mum  phase  reconstruction. 

With  an  adaptive  cepstral  gate  of  .8P  for  voiced  speech  and  a  fixed 
9  msec,  cepstral  gate  for  unvoiced  speech,  no  significant  differences  were 
noted  from  use  of  a  fixed  3.2  msec.  gate.  Consequently,  it  appears  thatany 
phase  errors  in  the  impulse  response  estimate  may  not  be  "noticed"  because 
of  the  overlap-addition  property  of  the  Hamming  window  discussed  In  section 
6.3.3,  and  since  the  windowed  waveform  has  a  low-time  complex  cepstrum. 

We  argued  In  Section  5.1.3  that  the  windowing  procedure  used  here  may 
not  be  suitable  for  spectral  magnitude  estimation.  From  Fig.  5.1  linear 
interpolation  of  complex  harmonic  samples  possibly  results  In  a  null  In  the 
magnitude  between  harmonics  when  the  phase  increment  Is  near  ir.  When  the 
pitch  estimate  Is  Inaccurate,  these  nulls  may  In  fact  be  sampled  upon  recon- 


structlon.  This  unfortunate  situation  suggests  that  the  magnitude  should 
perhaps  be  obtained  In  a  different  manner  from  that  of  the  phase.  The 
spectral  envelope  system  to  be  discussed.  In  fact,  applies  this  philosophy. 

In  spite  of  this  possibility,  however,  we  have  found  that  the  use  of 
a  Hanning  window,  of  duration  twice  a  pitch  period,  within  a  minimum  phase 
system  derived  from  the  magnitude  only,  generates  higher  quality  synthetic 
speech  than  a  system  where  a  longer  window  of  duration  40  msec* Is  applied 
[28].  A  longer  window  Introduces  other  forms  of  degradation  which  have  been 
Investigated  by  the  author  In  a  study  conducted  In  parallel  with  this  thesis 
[28]. 

6.4  Spectral  Envelope  Speech  Analysis-Synthesis 

In  this  section  we  first  briefly  review  the  structure  of  the  spectral 
envelope  speech  analysis-synthesis  system  proposed  by  Paul [26].  This  system 
relies  on  a  magnitude  estimate  only  and  Is  based  on  the  harmonic  samples  of 
the  Fourier  transform  of  the  windowed  speech  waveform.  A  phase  estimate  of 
the  vocal  tract  system  function  is  then  Introduced  by  the  linear  Interpola¬ 
tion  scheme  of  section  5.1.1,  also  based  on  harmonic  samples.  This  phase 
estimation  procedure  Is  therefore  easily  appended  to  Paul's  original  system. 

6,4.1  Minimum  Phase  Analysis-Synthesis 

Conceptually,  the  minimum  phase  spectral  envelope  analysis-synthesis 
system  Is  similar  to  the  homomorphic  configuration  of  section  6.3.1.  The 
primary  difference  lies  In  the  magnitude  estimation  procedure.  An  outline 
of  this  scheme  Is  given  as  follows: 

(I)  Apply  a  Hamming  window  to  the  speech  waveform. 

(II)  Find  the  location  and  value  of  the  harmonic  peaks  of  the  magnl- 


126 


tude  (or  logmagnltude)  of  the  resulting  Fourier  transform 
by  peak  picking.  A  pitch  adaptive  peak  picking  algorithm 
Is  used. 

(Ill)  Linearly  interpolate  the  harmonic  samples  of  the  magnitude 
(or  lo<pagn1tude). 

A  minimum  phase  estimate  Is  obtained  by  applying  the  causal  sequence 
of  the  form  In  (4.13)  to  the  real  cepstrum. 

6.4.2  Phase  Envelope  Estimation  By  Linear  Interpolation 

Step  (11)  above  gives  the  locations  of  harmonic  peaks.  We  can  use 
these  locations  to  find  the  values  of  H(w)  at  harmonics.  Given  these  values 
we  then  directly  apply  the  linear  Interpolation  scheme  of  section  5.1.1  to 
obtain  a  phase  estimate  which  preserves  the  desired  phase  envelope. 

We  must  of  course  align  the  waveform  so  that  condition  (5.16)  Is  sat¬ 
isfied.  In  particular,  the  heuristic  pre-alignment  process  of  section  6.3.3 
Is  utilized.  The  precise  duration  of  the  window,  however.  Is  not  crucial 
since  we  do  not  rely  on  the  window  itself  to  perform  Interpolation  In  the 
frequency  domain.  In  our  analysis,  we  choose  a  window  length  of  twice  the 
pitch  period  which  approximately  preserves  the  harmonic  values  and  Is  short 
enough  so  that  statlonarlty  of  the  speech  waveform  approximately  holds. 

6.4.3  Short-Time  Reconstruction 

With  our  analysis  scheme,  the  synthesizer  of  section  6.3.1  will  In 
theory  recover  a  periodic  waveform  exactly,  and  thus  acts  as  an  Identity 
system  for  such  waveforms.  This  Is  because  the  analysis  preserves  the 
harmonic  values  which  completely  characterizes  a  periodic  waveform.  Note 


127 


that  whan  pre-alignment  Is  such  that  condition  (5.16)  Is  not  satisfied,  but 
tha  alignment  Is  consistent  from  frame  to  frame,  the  Identity,  nevertheless, 
holds. 

6.4.4  Informal  Listening  Tests 

When  compared  with  Its  minimum  phase  counterpart,  the  system  with  mixed 
phase  provided  small,  but  audible.  Improvement  In  quality.  The  two  systems 
were  Identical  In  both  the  analysis  and  synthesis  stages  except  for  the 
phase  estimation  procedure  in  the  analysis  and  the  elimination  of  post¬ 
alignment  in  the  minimum  phase  version.  In  both  schemes  the  magnitude  esti¬ 
mate  Is  derived  by  linear  Interpolation  of  samples  of  the  magnitude  of  H(u) 
at  harmonics. 

The  informal  A/B  listening  test  of  section  6.3.4  was  performed.  The 
system  with  mixed  phase  was  judged  roughly  532  of  the  time  to  be  closer  In 
quality  to  the  original  than  Its  minimum  phase  counterpart,  the  minimum  phase 
version  102  of  the  time  was  judged  closer,  and  372  of  the  time  the  two  were 
Indistinguishable. 


CHAPTER  7 

SUMMARY  ANO  FUTURE  RESEARCH 


In  this  chapter  we  summarize  the  major  results  of  the  thesis  and  dis¬ 
cuss  a  number  of  directions  for  future  research. 

7,1  Summary 

In  this  dissertation  we  have  considered  techniques  of  phase  estimation 
for  the  purpose  of  speech  analysis  and  synthesis.  Our  methods  fall  rough¬ 
ly  within  the  categories  of  direct  and  Indirect  approaches. 

A  number  of  Indirect  techniques  generate  the  phase  of  a  Fourier  Trans¬ 
form  from  Its  magnitude  and  a  priori  knowledge  of  the  desired  sequence  or 
phase.  Both  closed  form  and  Iterative  solutions  were  developed  for  phase 
retrieval  from  magnitude  under  various  constraints.  An  alternate  Indirect 
means  of  estimating  phase  from  a  magnitude  relies  on  transformation  of  the 
speech  waveform  to  create  a  minimum  phase  Impulse  response  whose  magnitude 
Is  estimated.  An  Inverse  operation  provides  an  estimate  of  the  original 
phase  function. 

Direct  approaches  do  not  require  an  estimate  of  a  magnitude  function, 
but  require  samples  of  the  desired  system  function  or  partial  knowledge  of 
its  phase  derived  from  the  speech  waveform.  In  particular,  phase  estimates 
based  on  harmonic  samples  of  the  desired  system  function  were  incorporated 
within  two  speech  analysis-synthesis  systems  with  the  result  of  high-quality 
synthetic  speech. 

We  have  also  developed  an  Iterative  technique  for  magnitude  retrieval 
from  phase.  This  algorithm  provides  an  alternative  to  the  Hilbert  trans¬ 
form  for  obtaining  the  magnitude  from  the  phase  of  a  minimum  phase  sequence. 
The  technique  does  not  require  an  unwrapped  phase,  but  simply  Its  principal 
value.  The  Iteration  also  provides  a  means  for  magnitude  recovery  when 
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only  the  phase  of  a  finite  length  mixed  phase  sequence  Is  given.  Further¬ 
more,  the  Iterative  minimum  phase  signal  reconstruction  serves  as  the  major 
component  within  a  new  phase  unwrapping  algorithm  which  does  not  require 
modulo  2 it  considerations. 

This  study  represents  only  an  Initial  step  In  developing  techniques 
of  phase  estimation  In  the  realm  of  speech  analysis  and  perhaps,  equally 
Important,  In  other ‘areas  such  as  Image  restoration  where  an  accurate  phase 
estimate  may  be  particularly  significant. 

Below  we  outline  a  number  of  possible  directions  for  future  research. 

7.2  Suggestions  For  Future  Research 
Sensitivity  Analysis 

The  methods  of  phase  and  magnitude  retrieval  In  sections  3.2  and  4.2 
rely  on  exact  knowledge  of  a  magnitude  or  phase  function,  or  values  of  a 
sequence.  Sensitivity  of  either  the  closed  form  or  Iterative  solutions  to 
noise  or  other  degradation  Is  not  completely  understood.  A  number  of  obser¬ 
vations,  however,  can  be  made.  When  causality  and  phase  continuity  con¬ 
straints  are  Imposed,  clearly  the  uniqueness  argument  of  section  4.2  still 
holds.  That  Is,  given  a  degraded  phase  there  exists  one  degraded  magnitude 
which  corresponds  to  a  minimum  phase  sequence.  With  a  similar  argument,  we 
conclude  that  under  a  finite  length  constraint  there  exists  one  magnitude 
function  for  a  given  distorted  phase. 

In  fact,  with  phase  distortion,  our  Iterative  algorithm  to  recover 
magnitude  from  phase  converges  to  the  unique  predicted  magnitude.  The 
magnitude  function,  however,  may  differ  significantly  from  the  original 
distortionless  magnitude.  For  example.  In  minimum  phase  reconstruction 
by  Iteration,  a  maximum  phase  zero  lying  near  the  unit  circle  may  become 
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a  minimum  phase  zero  In  the  presence  of  noise  rather  than  a  minimum  phase 
pole  as  determined  by  (2.28).  This  phenomenon,  however.  Is  not  a  function 
of  the  iterative  procedure,  but  of  the  constraints  Imposed. 

The  sensitivity  of  the  dual  problem  of  recovery  of  phase  from  magni¬ 
tude  Is  more  difficult  to  characterize.  An  ambiguity  arises  since  addition 
of  noise  may  cause  the  breakdown  of  our  assumed  rational  models  and  con¬ 
ditions  of  section  3.2. 

Another  sensitivity  Issue  involves  the  effect  of  noise  or  other  degrada¬ 
tion  on  the  technique  of  section  4.1.4  where  a  phase  estimate  is  obtained 
by  transformation  of  a  quasi -period 1c  waveform.  We  Illustrated  one  example 
of  this  approach  to  Improve  the  phase  estimate  over  a  direct  procedure  In 
the  presence  of  noise.  However,  the  mapping  of  magnitude  degradation  through 
an  inverse  transformation  to  phase  degradation  Is  not  understood. 

Convergence  Issues 

The  iterative  algorithm  for  phase  and  magnitude  retrieval  were  found 
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in  practlve  to  generate  converging  solutions  sometimes  slowly  (e.g. ,  after 
several  hundred  iterations)  and  sometimes  quickly  (e.g.,  after  a  few  iter¬ 
ations).  Consequently,  determining  rates  of  convergence  In  terms  of  spec¬ 
tral  structure,  and  methods  for  quickening  convergence  are  useful  areas  of 
research. 

Another  question  Involves  the  existence  of  rigorous  proofs  demonstrat¬ 
ing  unique  convergence  of  015(a))  and  MjJid)  derived  from  the  Iterative  algor¬ 
ithms. 

Alternative  Constraints  for  Phase  Retrieval 


Given  a  magnitude  function,  we  might  consider  constraints  for  unambig¬ 
uous  phase  retrieval  where  specific  nonzero  values  of  h(n)  are  not  required. 
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Consider,  for  example,  the  minimum  and  maximum  values  of  h(n)  or  e^U). 

A  question  which  naturally  follows  is  whether  such  constraints  can  be 
imposed  within  an  Iterative  scheme. 

Partial  Magnitude  -  Partial  Phase 

In  sections  3.1  and  4.2.1  samples  of  the  magnitude  or  phase  along  with 
certain  nonzero  values  of  h(n)  were  sufficient  for  phase  or  magnitude  re¬ 
trieval.  One  obvious  extension  of  this  result  concerns  signal  reconstruc¬ 
tion  from  samples  of  both  magnitude  and  phase.  We  might  consider  possible 
ways  to  distribute  such  samples  to  unambiguously  recover  the  entire  magni¬ 
tude  and  phase  functions.  In  the  context  of  coding,  certain  sampling  dis¬ 
tributions  may  cause  little  degradation  of  the  reconstructed  signal. 

"Optimal11  Transformations 

In  section  4.3  we  have  only  touched  upon  the  use  of  transformations  in 
phase  estimation.  What  Is  needed  is  a  method  for  designing  "optimal" 
transformations  which  create  minimum  phase  sequences  with  specific  desir¬ 
able  properties.  For  example,  in  the  context  of  homomorphic  deconvolution 
the  real  cepstrum  of  the  modified  sequence  should  be  as  low  time,  as  pos¬ 
sible,  but  not  susceptible  to  quantization  noise  as  occurs  with  exponential 
weighting.  These  properties  are  also  useful  in  coding  samples  of  the  log- 
magnitude  function.  In  the  context  of  signal  enhancement,  the  transforma¬ 
tion  should  be  such  that  the  resulting  magnitude  function  Is  suitable  to 
estimation  by  a  specific  noise  reduction  scheme  such  as  spectral  subtraction. 

Phase  Modeling 

Suppose  an  accurate  estimate  of  the  unwrapped  phase  Is  available.  If 
this  phase  as  well  as  magnitude  function  is  to  be  efficiently  transmitted 
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in,  for  example,  a  speech  communications  system,  a  small  parameter  set 
must  be  derived  which  represents  the  entire  phase  function.  When  the  un¬ 
wrapped  phase  Is  bandllmited,  samples  of  the  phase  are  adequate.  Alter¬ 
nately,  we  may  consider  ways  to  directly  model  the  phase,  or  phase  deriva¬ 
tive,  by  a  rational  function. 

A  Phase-Only  Vocoder 

We  saw  In  section  6.3.3  that  the  homomorphic  analysis-synthesis  scheme 
is  an  identity  system  with  respect  to  a  periodic  waveform  when  an  adaptive 
Hamming  window  Is  applied.  The  quality  of  the  overall  system  does  not  sig¬ 
nificantly  differ  when  a  low-time  or  unity  cepstral  gate  is  applied  since 
the  complex  logarithm  is  nearly  bandlimited.  Aligning  the  Hamming  window 
and  adapting  Its  duration  to  the  pitch  period  yields  logmagnitude  and  un¬ 
wrapped  phase  functions  whose  slowly  varying  components  approximately  re¬ 
present  the  windowed  waveform. 

Since  the  windowed  waveform  is  of  finite  extent,  from  section  4.2  the 
phase  Is  sufficient  for  complete  waveform  characterization.  Therefore, 
coding  samples  or  a  small  set  of  model  parameters  of  the  bandlimited  phase 
function  may  be  adequate  for  signal  reconstruction. 
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